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Abstract 

We  study  the  dynamics  and  properties  of  a  turbulent  flame,  formed  in  the  presence  of  subsonic,  high-speed,  homoge¬ 
neous,  isotropic  Kolmogorov-type  turbulence  in  an  unconfined  system.  Direct  numerical  simulations  are  performed 
with  Athena-RFX,  a  massively  parallel,  fully  compressible,  high-order,  dimensionally  unsplit,  reactive-flow  code.  A 
simplified  reaction-diffusion  model  represents  a  stoichiometric  H2-air  mixture.  The  system  being  modeled  represents 
turbulent  combustion  in  the  well-stirred  reactor  regime,  with  Damkohler  number  Da  =  0.1  and  the  turbulent  velocity 
at  the  energy  injection  scale  30  times  larger  than  the  laminar  flame  speed.  The  simulations  show  that  flame  interaction 
with  high-speed  turbulence  forms  a  steadily  propagating  turbulent  flame  with  a  flame  width  approximately  twice  the 
energy  injection  scale  and  a  speed  four  times  the  laminar  flame  speed.  A  method  for  reconstructing  the  internal  flame 
structure  is  described  and  used  to  show  that  the  turbulent  flame  consists  of  tightly  folded  flamelets.  The  internal  struc¬ 
ture  of  these  is  virtually  identical  to  that  of  the  planar  laminar  flame  with  the  preheat  zone  broadened  by  approximately 
a  factor  of  two.  The  turbulent  cascade  fails  to  penetrate  the  internal  flame  structure,  and  so  the  action  of  small-scale 
turbulence  is  suppressed  throughout  most  of  the  flame.  Finally,  our  results  suggest  that  for  stoichiometric  H2-air 
mixtures  any  substantial  flame  broadening  by  the  action  of  turbulence  cannot  be  expected  in  all  subsonic  regimes. 

Key  words:  Turbulent  premixed  combustion.  Turbulence,  Flamelet,  Flame  structure.  Hydrogen 


1.  Introduction 

In  the  last  twenty  years,  studies  of  premixed  turbulent  combustion  have  witnessed  an  explosive  growth.  The 
motivation  for  this  comes  from  a  remarkably  broad  range  of  both  engineering  and  basic  science  applications  -  from 
the  design  of  internal  combustion  engines,  to  problems  of  ensuring  safety  in  mines,  to  the  dynamics  and  properties  of 
turbulent  thermonuclear  flames  powering  type  la  supernovae.  Such  rapid  development  has  been  enabled  by  significant 
advances  in  experimental  techniques,  in  particular  in  terms  of  the  diagnostics,  as  well  as  by  the  substantial  increase 
in  the  capabilities  of  the  computational  infrastructure  and  numerical  algorithms,  which  made  the  direct  numerical 
simulation  (DNS)  a  viable,  and  often  indispensable,  tool  for  turbulent  combustion  research. 

Despite  the  significant  increase  in  the  body  of  experimental  and  numerical  data,  the  overall  paradigm  still  prevail¬ 
ing  today  is  reflected  in  a  variety  of  combustion  regime  diagrams  [1,2,  3],  which  attempt  to  provide  a  comprehensive 
classification  of  the  different  modes  of  turbulence-flame  interaction.  The  key  underlying  assumption  is  that  such  clas¬ 
sification  can  be  made  based  on  a  very  limited  set  of  parameters,  namely,  characteristic  large-scale  turbulent  velocity 
and  spatial  scale.  With  this  ansatz,  the  diagrams  can  be  constructed  by  comparing  various  turbulent  timescales,  e.g., 
eddy  turnover  timescales  on  the  integral  and  Kolmogorov  scales,  with  those  of  the  unperturbed  laminar  flame.  The 
resulting  classifications  include  such  regimes  as  “wrinkled”  and  “corrugated  flamelets,”  “distributed  reaction  zones,” 
and  “well-stirred  reactors”,  which  suggest  the  typical  qualitative  structure  of  the  turbulent  flame  under  specific  condi¬ 
tions. 

The  basis  for  this  paradigm  is  the  picture  which  was  first  suggested  by  Damkohler  almost  70  years  ago  [4,  5,  6]. 
Farge-scale  motions  are  responsible  for  the  overall  stretching  and  folding  of  the  flame  that  constitute  the  flame  brush. 
This  process  increases  the  flame  surface  area,  which  directly  determines  the  global  properties  of  the  turbulent  flame, 
such  as  its  width  and  speed.  Therefore,  depending  on  whether  the  characteristic  timescale,  associated  with  those 
large-scale  motions,  is  greater  or  smaller  than  the  laminar  flame  self-crossing  time,  the  flame  is  either  able  or  not  to 
reorganize  and  adjust  itself  to  the  action  of  turbulence. 

On  the  other  hand,  small  scales  penetrate  inside  individual  flamelets,  thus  affecting  their  internal  structure  and 
broadening  their  preheat  and  reaction  zones.  The  efficiency  of  this  process  is  controlled  by  the  magnitude  of  the 
Kolmogorov  scale  with  respect  to  the  laminar  flame  width,  6l,  or  the  size  of  the  reaction  zone  in  the  laminar  flame. 
Once  the  Kolmogorov  scale  becomes  smaller,  it  is  hypothesized  that  the  turbulent  cascade  is  able  to  penetrate  the 
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internal  structure  of  the  flamelets.  As  a  result,  turbulent  transport  becomes  comparable  to,  or  exceeds,  molecular 
diffusive  processes,  and  the  flamelet  width  and  burning  velocity  increase. 

Overall,  in  this  picture  the  actions  of  large  and  small  scales  are  quite  distinct.  The  large  scales  determine  the  global 
turbulent  burning  speed,  which  increases  compared  to  the  laminar  burning  speed  in  proportion  to  the  increase  in  the 
flamelet  surface  area.  At  the  same  time,  small  scales  affect  the  turbulent  flame  speed  only  by  modifying  the  local 
burning  velocity  of  individual  flamelets. 

This  dichotomy  in  the  action  of  the  large  and  small  scales  by  and  large  has  its  root  in  the  behavior  of  a  passive 
scalar  in  turbulent  nonreactive  flows.  The  fact  that  in  non-reacting  turbulence  the  inertial  range  extends  to  scales 
smaller  than  6l,  however,  does  not  automatically  mean  that  turbulent  motions  of  the  same  intensity  penetrate  inside 
the  flamelets,  nor  does  it  mean  that  they  have  the  same  effect  on  the  internal  flamelet  structure  as  they  would  have  on 
a  passively  advected  scalar. 

The  dominant  majority  of  research  efforts  up-to-date  has  focused  on  regimes  characterized  by  smaller  ratios  of 
turbulent  velocities  to  the  laminar  flame  speeds  (see  reviews  by  [5,  6]),  where  there  are  the  most  applications  of 
practical  interest.  In  these  cases,  penetration  by  small-scale  turbulent  motions  into  the  internal  flamelet  structure  is 
either  completely  suppressed  ,or  it  is  heavily  dominated  by  flame  wrinkling  due  to  large-scale  motions.  Results  of 
these  efforts  are  in  agreement  with  the  Damkohler  concept  discussed  above,  and  they  support  the  classification  of  such 
regimes  as  “wrinkled”  or  “corrugated”  flamelets. 

The  investigation  of  the  regimes  in  which  turbulent  velocities  are  significantly  larger  than  the  laminar  flame  speed, 
S  L,  on  all  scales,  including  that  of  the  laminar  flame  thickness,  presents  a  number  of  both  experimental  and  numerical 
challenges.  Hereafter,  we  refer  to  this  mode  as  high-speed  turbulent  combustion,  which  encompasses  regimes  such 
as  thickened  reaction  zones,  distributed  flames,  and  well-stirred  reactors.  These  are  the  regimes  in  which  substantial 
flame  broadening  by  turbulent  transport  has  been  hypothesized.  We  also  assume  that  turbulence  is  the  only  process 
that  can  affect  the  structure  and  properties  of  the  flame.  We  do  not  consider  situations  in  which  the  flame  is  altered  or 
broadened  by  fuel  preconditioning,  compression  of  the  overall  system,  or  propagation  of  large-scale  shocks. 

Probing  such  regimes  experimentally  requires  either  high  turbulent  intensities  or  low  laminar  flame  speeds.  Gen¬ 
erating  and  sustaining  high-speed  reactive  turbulence  with  desired  properties  may  be  difficult  in  a  laboratory  setting. 
More  important,  however,  is  that  if  5  z,  is  high,  the  maximum  ratios  of  the  integral  turbulent  velocity  Uito  S  l  cannot 
be  too  large,  or  the  overall  flow  will  no  longer  be  subsonic.  Therefore,  almost  all  studies  that  were  able  to  achieve 
Ui/S L  ^  10  -  20,  used  lean  hydrocarbon  fuels,  which  have  low  values  of  Sl  [7,  8,  9]  (also  see  the  review  by  [5]). 

Of  particular  interest  is  the  work  by  Dunn  et  al.  [10],  in  which  a  piloted  premixed  jet  burner  was  used  to  achieve 
UilS L  in  the  range  40  -  390,  corresponding  to  Karlovitz  numbers  Ka  =  100  -  3500  and  Damkohler  numbers  Da  - 
0.069  -  0.0053.  In  their  lowest-velocity  case,  which  would  be  characterized  as  being  well  into  the  distributed-flame 
regime,  there  was  no  evidence  of  flame  broadening  based  on  the  temperature  and  OH  images.  This  is  in  contradiction 
to  predictions  of  the  traditional  classification  of  the  combustion  regimes.  Substantially  higher  turbulent  intensities 
were  required  to  produce  a  flow  structure,  which  was  characterized  by  the  authors  as  a  distributed  mode  of  burning, 
although  this  determination  was  done  mostly  on  the  qualitative  grounds. 

In  numerical  simulations,  generating  turbulence  of  arbitrary  intensity  does  not  present  any  complications.  Instead, 
modeling  high-speed  turbulent  reactive  flows  faces  different  challenges  which  can  be  as  restrictive  as  those  encoun¬ 
tered  in  experiments.  One  major  difficulty  is  that  the  Kolmogorov  scale  can  be  substantially  smaller  than  6l,  thereby 
making  fully-resolved  DNS  difficult  even  with  the  largest  computational  resources.  Furthermore,  when  fluid  velocities 
are  high  compared  with  the  speed  of  sound,  it  is  necessary  to  use  fully  compressible  integration  algorithms,  for  which 
the  numerical  time  step  is  limited  by  the  CFL  condition  further  increasing  the  computational  cost. 

Several  studies  of  chemical  flames  achieved  values  of  UifSi  in  the  range  3  -  10  [11,  12,  13]  in  a  flat-flame 
configuration  and  UijS l  -  20  in  the  case  of  a  bunsen  flame  [14]  (also  see  the  review  by  [5]).  In  all  of  these  studies, 
the  reaction  zone  has  a  thin  sheet-like  structure.  Moreover,  Bell  et  al.  [12]  demonstrate  that  the  increase  in  turbulent 
velocity  can  almost  completely  be  accounted  for  by  the  corresponding  increase  in  the  flamelet  surface  area.  This  is 
the  case  even  though  the  turbulent  flame  speed  appears  to  be  growing  somewhat  faster  than  the  surface  area.  It  should 
be  noted,  that  studies  considering  the  flat-flame  configuration  [11,  12,  13]  used  decaying  turbulence,  which  precluded 
the  system  from  reaching  a  steady  state. 

Aspden  et  al.  [15]  reported  flame  broadening  and  formation  of  a  distributed  reaction  zone  in  a  DNS  under  the 
action  of  fast  turbulence.  This  study  considers  the  interaction  of  a  thermonuclear  flame,  driven  by  reactions, 

with  steady  driven  turbulence  under  the  conditions  characteristic  of  the  white  dwarf  interior  during  the  late  stages  of 
the  type  la  supernova  explosion.  In  the  fastest  turbulence  case,  the  value  of  UijS l  ~  70  was  achieved  corresponding 
to  Ka  -  228.  In  this  simulation,  formation  of  a  broad  turbulent  flame  was  observed  with  complex  temperature  and 
burning-rate  structure.  Moreover,  the  joint  PDFs  of  various  quantities,  such  as  carbon  mass  fraction  and  temperature, 
were  characteristic  of  a  Lewis  number  Le  =  1  system,  whereas  the  actual  values  of  Le  are  extremely  large.  This  fact, 
along  with  wide  smooth  averaged  profiles  of  various  quantities  throughout  the  turbulent  flame,  led  to  the  conclusion 
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that  the  resulting  flame  is  dominated  by  turbulent  transport  and  presents  an  example  of  a  broad  distributed  reaction 
zone. 

Even  given  the  efforts  discussed  above,  the  physics  of  high-speed  turbulent  combustion  still  has  many  unanswered 
questions.  With  the  exception  of  the  experimental  results  of  [10]  and  numerical  results  of  [15],  evidence  of  significant 
flame  broadening  and  the  formation  of  distributed  reaction  zones  is  scarce  [5].  Moreover,  results  of  all  studies  suggest 
that  such  broadening,  even  if  does  occur,  requires  substantially  higher  turbulent  intensities  than  traditionally  predicted. 
Whether  the  effects  of  turbulence  are  indeed  so  much  weaker  than  can  be  expected  based  on  simple  physical  grounds, 
as  well  as  the  mechanisms  suppressing  the  action  of  turbulence,  remain  unclear. 

The  work  by  [10]  and  [15]  emphasizes  the  importance  of  the  detailed  and  reliable  diagnostics.  While  both  studies 
present  interesting  experimental  and  numerical  investigations  of  the  previously  unexplored  regimes,  they  both  rely 
on  indirect  evidence  in  their  determination  of  the  observed  flame  broadening.  Of  paramount  importance  are  the 
measurements  that  would,  for  example,  directly  answer  questions  regarding  the  actual  multidimensional  structure  of 
the  region  in  which  most  of  the  energy  release  takes  place,  the  nature  of  flame  broadening,  and  the  effects  of  turbulence 
on  different  regions  of  the  flame. 

Our  goal  in  this  work  is  to  address  several  key  questions  concerning  the  dynamics  and  properties  of  the  turbulent 
flames  formed  under  the  action  of  high-speed  turbulence.  First  of  all,  it  is  important  to  understand  (1)  the  global 
properties  of  the  turbulent  flame  brush  in  this  regime,  i.e.,  its  width,  speed,  etc.;  (2)  the  internal  structure  of  the  flame 
brush;  and  (3)  the  internal  structure  of  the  flamelets  folded  inside  the  flame  brush,  if  such  flamelets  can  be  identified  at 
all.  It  is  important  to  address  the  last  two  points  quantitatively  through  direct  measurement.  All  of  these  characteristics 
are  determined  by  turbulence-flame  interactions,  and  the  resulting  system  is  the  product  of  the  collective  action  of  the 
full  turbulent  energy  cascade.  Therefore,  it  is  also  very  important  to  determine  the  relative  role  of  large-  and  small- 
scale  motions  on  the  global  and  local  properties  of  the  turbulent  flame. 

In  this  paper,  we  begin  to  address  these  questions  by  presenting  three  simulations,  which  are  designed  to  represent 
the  turbulent  combustion  process  in  an  unconfined  (ideally  infinite)  space.  Initially,  the  flame  is  a  planar  surface 
separating  half  of  the  domain  containing  fuel  from  the  half  of  the  domain  containing  products  and  immersed  in  a 
high-intensity  turbulent  flow  field  with  the  Kolmogorov-type  spectrum  [16].  It  is  assumed  that  the  energy  injection 
scale,  as  well  as  the  turbulent  integral  scale,  are  both  finite  and  much  smaller  than  the  overall  system  size.  The 
evolution  of  such  a  system  represents  the  evolution  of  an  infinite  globally  planar  turbulent  flame  brush.  Such  an 
idealized  setting  allows  us  to  exclude  effects  such  as  walls,  boundaries,  and  the  system  geometry  from  consideration 
and  to  focus  exclusively  on  the  role  played  by  turbulence. 

We  consider  driven  turbulence  of  sufficiently  high  intensity  to  place  the  system  well  into  the  well-stirred  reactor 
regime  according  to  the  traditional  classification.  We  consider  the  stoichiometric  reactive  mixture  to  minimize  the 
thermo-diffusive  effects  as  well  as  the  possibility  of  flame  extinction  under  the  action  of  intense  turbulence.  We 
use  a  simplified  reaction  model  designed  to  accurately  represent  the  laminar  flame  and  detonation  properties  of  the 
stoichiometric  H2-air  mixture.  A  high  laminar  flame  speed  of  the  fuel  leads  to  turbulent  velocities  which  are  relatively 
large  compared  with  the  sound  speed  in  cold  fuel.  The  resulting  flow  can  no  longer  be  considered  incompressible, 
which  requires  the  use  of  a  high-order,  dimensionally  unsplit,  fully  compressible  integration  method. 

The  three  simulations  on  which  we  focus  in  this  paper  are  part  of  a  larger  series  of  numerical  models.  Their  goal 
is  to  survey  the  full  range  of  subsonic  high-speed  turbulent  combustion  regimes  in  a  variety  of  the  reactive  mixtures. 
This  work  is  intended  to  identify  the  overall  framework  for  analysis  of  such  numerical  models. 


2.  Formulation  of  the  Problem 


2.7.  Physical  model 

We  solve  the  system  of  unsteady,  compressible,  reactive  flow  equations, 

dp 

|.V.(pu)  = 


dt 


-I-  V  ■  (pu  ®  u)  H-  VP  = 
ot 

-H  V  ■  ((£  -H  P)u)  -H  V  ■  (TfVP)  = 
^  -H  V  ■  (pTu)  -I-  V  ■  (pDVy)  = 


0, 

(1) 

0, 

(2) 

pqw. 

(3) 

pw. 

(4) 

Here  p  is  the  mass  density,  u  is  the  velocity,  E  is  the  energy  density,  P  is  the  pressure,  Y  is  the  mass  fraction  of  the 
reactant,  K  and  D  are  the  coefficients  of  thermal  conduction  and  molecular  diffusion,  q  is  the  chemical  energy  release. 
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and  w  is  the  reaction  source  term.  The  equation  of  state  is  that  of  an  ideal  gas.  We  do  not  include  physical  viscosity  in 
our  model.  Therefore,  dissipation  is  provided  by  the  numerical  viscosity.  This  issue,  and  the  way  we  use  the  inherent 
numerical  viscosity  in  the  algorithm,  will  be  discussed  in  detail  in  §  2.4.  Finally,  since  our  method  of  turbulence 
driving,  described  in  §  2.2,  does  not  involve  modifying  the  fluid  equations,  eqs.  (1)  -  (4)  do  not  contain  any  forcing 
term. 

Chemical  reactions  are  modeled  using  the  first-order  Arrhenius  kinetics 

dY  I  Q\ 

-  =  »  =  -Apfexp(--).  (5) 

where  A  is  the  pre-exponential  factor,  Q  is  the  activation  energy,  and  R  is  the  universal  gas  constant. 

We  assume  that  both  thermal  conduction  and  species  diffusion  have  similar  temperature  dependence 

'j^n  ^  Y'n 

D  =  Dq  — ,  =  kq  — ,  (6) 

P  pCp  p 

where  Do,  /co,  and  n  are  constants,  and  Cp  -  yRjMij  -  1)  is  the  specific  heat  at  constant  pressure.  The  Lewis  number 
Le  -  /Co /Do  is  set  to  be  equal  to  one  and  it  is  independent  of  the  thermodynamic  conditions  of  the  flow.  External 
forces,  the  Soret  and  Dufour  effects,  pressure  gradient  diffusion  as  well  as  the  effects  of  radiation  are  assumed  to  be 
negligible.  The  last  assumption,  however,  requires  further  verification.  Contributions  of  the  molecular  dissociation, 
atomic  ionization,  and  radiative  heat  losses  to  the  equation  of  state  are  assumed  to  be  absorbed  by  the  model  value  of 
the  adiabatic  index  y. 

Reaction  model  parameters  are  summarized  in  table  1.  They  are  based  on  the  simplified  reaction  model  of  [17] 
designed  to  represent  the  stoichiometric  H2-air  mixture.  The  values  of  the  transport  coefficients  are  comparable  to 
those  of  air  at  the  same  conditions.  This  reaction  model  reproduces  the  key  characteristics  both  of  the  laminar  flames 
and  the  multidimensional  detonations  in  the  given  reactant  mixture,  such  as  the  laminar  flame  width,  speed,  detona¬ 
tion  velocity  and  detonation  cell  size,  etc.  as  well  as  the  dependence  of  these  quantities  on  pressure  and  temperature. 
It  has  also  been  demonstrated  to  provide  good  qualitative  agreement  with  experimental  data  in  more  complex  ap¬ 
plications  involving  flame  acceleration  and  deflagration-to-detonation  transition  in  channels  with  obstacles  [17,  18]. 
These  properties  of  the  reaction  model,  along  with  its  low  computational  cost,  make  it  a  practical  choice  for  large- 
scale  multidimensional  simulations  aimed  at  qualitative  analysis  of  the  turbulence-flame  interaction  in  stoichiometric 
hydrogen-air  mixtures  and  other  similar  fuels.  It  is  important  to  keep  in  mind  that  the  first-order  Arrhenius  kinetics 
is  not  able  to  capture  the  full  complexity  of  hydrogen  combustion.  Therefore,  certain  care  must  be  exercised  when 
comparing  results  obtained  using  this  model  with  the  actual  experimental  data. 

2.2.  Numerical  method 

To  model  the  process  of  turbulence-flame  interaction,  we  use  the  code  Athena-RFX  [19]  -  the  reactive  flow 
extension  of  the  magnetohydrodynamic  code  Athena  [20,  21,  22].  Athena-RFX  is  a  fixed-grid  massively  parallel 
code.  It  implements  several  higher-order  fully  conservative  Godunov-type  methods  for  integration  of  fluid  equations. 
In  this  work  we  use  the  method  based  on  the  fully  unsplit  comer  transport  upwind  (CTU)  algorithm  of  [23]  and  its 
3D  extension  presented  by  [24].  In  particular,  Athena  implements  a  variant  of  this  method  [21],  which  requires  six 
Riemann  solves  per  cell  instead  of  the  twelve  in  the  original  method  of  [24].  This  integration  scheme  uses  PPM 
spatial  reconstmction  [25]  in  conjunction  with  the  approximate  nonlinear  HLLC  Riemann  solver  to  achieve  3rd- 
order  accuracy  in  space  and  2nd-order  accuracy  in  time.  A  detailed  description  and  an  extensive  suite  of  tests  of  the 
hydrodynamic  integration  algorithm  and  its  implementation  in  Athena  can  be  found  in  [20,  21]. 

Diffusive  transport  is  incorporated  into  the  hydrodynamic  integration  algorithm  by  calculating  net  diffusive  fluxes 
across  each  cell  along  each  dimension.  This  is  done  using  second-order  finite  differences  with  flux  matching  to  ensure 
that  the  conservativity  is  preserved.  Those  fluxes  are  then  added  to  the  hydrodynamic  fluxes  prior  to  the  final  update 
of  the  state  vector  in  each  cell.  Source  terms,  describing  chemical  reactions,  are  coupled  using  Strang  splitting  to 
ensure  the  second-order  accuracy  of  the  overall  solution. 

The  overall  resulting  solver  is  formally  second-order  accurate  and  it  is  capable  of  providing  the  accuracy  of 
the  planar  laminar  flame  solution  of  <  1  -  2%  even  with  the  resolution  of  ~  4  cells  per  6l-  We  refer  to  [19]  for 
further  details  of  the  implementation  of  the  reactive-diffusive  extensions  in  Athena-RFX,  as  well  as  the  results  of  tests 
including  convergence  studies. 

2.3.  Turbulence  driving  method 

We  investigate  the  process  of  flame  interaction  with  steady  homogeneous  isotropic  turbulence,  described  by  the 
classical  Kolmogorov  theory  [16]  (hereafter  K41).  In  the  absence  of  persistent  energy  injection  into  the  system  at  a 
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large  scale,  turbulence  decays  on  a  characteristic  time  scale  of  the  large-scale  eddy  turnover  time.  Under  conditions 
considered  later  in  this  paper,  this  time  scale  is  almost  four  times  shorter  than  the  laminar  flame  self-crossing  time 
^lIS L  (see  tables  1  and  3).  Consequently,  in  order  to  study  the  quasi-steady  regime  of  turbulence-flame  interaction, 
the  flow  must  be  constantly  stirred  at  the  largest  scale  in  order  to  create  and  maintain  a  steady  energy  cascade  to 
smaller  scales. 

We  use  a  spectral  turbulence-driving  method  similar  to  the  one  used  in  [26,  27].  For  completeness  here  we 
summarize  the  key  algorithmic  stages  of  the  method.  More  detailed  discussion  of  the  method,  along  with  the  analysis 
of  its  properties  and  an  extensive  suite  of  tests  demonstrating  its  performance,  can  be  found  in  [19]. 

As  the  first  step,  we  consider  Fourier  transforms  of  the  velocity  perturbations  (5u'(k),  with  each  component  6u'- 
being  an  independent  realization  of  a  Gaussian  random  field  with  zero  mean  and  unit  variance.  Then  a  given  isotropic 
energy  injection  spectrum  S&(k)  is  superimposed  on  i5u'(k) 


6ui(k)  =  ^-j^6u'(k),  (7) 

where  k  -  |k|.  In  principle,  the  amplitude  and  phase  of  each  component  Sui  at  every  point  k  can  be  independently 
adjusted  to  produce  an  energy-injection  pattern  of  arbitrary  complexity. 

In  the  computational  domain  of  volume  V  =  L1L2T3,  where  Li  <  L2  <  L3,  discrete  Fourier  space  spans  wavenum¬ 
bers  kilAki  e  [-fN,fN),  where  i  6  1,2,3,  Ak,-  =  InlLi,  Ax  is  the  cell  size,  and  fy  -  LijlAx  -  Nil2  is  the  Nyquist 
frequency.  In  this  work  we  inject  energy  only  at  the  scale  L  -  Li  to  obtain  the  K41-type  spectrum.  In  particular. 


6S(k)  = 


0  V  k  =  jki,k2,k3)  ;  =  ^2  =  ^3  =  0,  ^  0  or  1, 

1  V  k  =  [ki,k2,k3):  ^(^)  =  0  or  1. 


(8) 


Finally,  non-solenoidal  components  of  the  velocity  perturbations  are  removed  using  the  orthogonal  projection  opera¬ 
tor,  which  ensures  that  V  ■  bu(x)  =  0. 

An  inverse  Fourier  transform  of  <5u(k)  gives  bu(x),  the  velocity  perturbation  field  in  the  physical  space.  Then  the 
bu(x)  are  added  to  the  velocity  field  u(x)  in  the  domain  on  every  time  step,  after  first  being  normalized  to  ensure 
constant  rate  e  of  kinetic-energy  injection.  Moreover,  the  total  momentum  in  the  perturbation  field  is  subtracted  from 
bu(x)  to  ensure  that  no  net  momentum  is  added  to  the  domain,  i.e.,  J  p6u  -  0.  The  overall  perturbation  pattern  is 
regenerated  at  every  time  interval  Afi,^  ss  5Ax/cj,  where  Cj  is  the  sound  speed  in  the  domain. 

This  turbulence  driving  method  produces  statistically  steady  turbulent  flows  with  energy  spectra  of  arbitrary  com¬ 
plexity  [19].  In  particular,  it  is  possible  to  obtain  K41-type  turbulence  with  the  inertial  range  of  the  energy  cascade 
extending  all  the  way  to  the  energy  injection  scale.  This  is  particularly  important,  given  the  limited  range  of  spatial 
scales  typically  accessible  in  a  simulation.  The  method  does  not  introduce  any  large-scale  anisotropies  or  global  flows. 
Moreover,  since  the  velocity  perturbation  field  is  purely  solenoidal,  no  compressions  or  rarefactions  are  artificially  in¬ 
duced  as  a  result  of  driving.  This  is  particularly  important  in  the  case  of  the  reactive  flows,  in  which  the  rate  of  energy 
generation  can  be  very  sensitive  to  the  local  variations  in  temperature  and  pressure. 


2.4.  Problem  setup  and  summary  performed  simulations 

Table  2  summarizes  key  parameters  of  the  performed  simulations.  The  main  difference  between  the  three  models 
is  the  resolution,  which  progressively  increases  from  Ax  =  in  SI  to  Ax  =  6lI32  in  S3.  One  of  the  principal  goals 
of  this  work,  discussed  in  §  1,  is  to  investigate  the  effects  of  small-scale  motions  and,  thereby,  to  differentiate  them 
from  the  role  played  by  large  scales.  The  most  direct  way  to  achieve  this  is  by  varying  the  viscosity.  In  particular, 
changing  the  resolution  changes  the  numerical  dissipation,  and  this  causes  the  effective  viscosity  in  the  domain  to 
vary.  This  allows  us  to  control  the  amount  of  kinetic  energy  contained  on  scales  comparable  to  or  smaller  than  6l, 
while  preserving  the  energy  spectrum  on  larger  scales. 

Figure  1  shows  instantaneous  spectra  in  the  domain  prior  to  ignition  in  simulations  SI  -  S3.  These  spectra 
represent  the  flow  field  in  the  nonreactive  turbulence.  While  the  spectra  are  virtually  the  same  on  scales  >  26l,  they 
vary  greatly  in  energy  contained  on  small  scales.  In  the  least  resolved  case,  SI,  motions  on  all  scales  <  6l  are  affected 
by  dissipation.  In  the  most  resolved  case,  S3,  the  inertial  range  extends  to  scales  which,  in  the  reactive  flow,  would 
be  deep  inside  the  flame.  As  a  result,  energy  on  the  small  scales  varies  by  up  to  a  factor  of  30  between  these  two 
cases.  By  comparing  the  flow  structure  in  these  simulations,  it  is  then  possible  to  establish  the  effects  of  small-scale 
turbulence  on  the  turbulent  flame.  In  particular,  this  allows  us  to  determine  the  degree  to  which  turbulent  cascade  is 
able  to  penetrate  and  alter  the  internal  structure  of  the  flame. 

A  similar  analysis  could  be  performed  by  including  physical  viscosity  in  the  model  and  by  varying  its  strength 
parametrically.  Such  an  analysis,  however,  would  not  represent  the  behavior  of  an  actual  fuel  any  more  than  the  one 
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which  relies  on  numerical  viscosity.  Moreover,  there  are  several  practical  complications  associated  with  using  physical 
viscosity  instead  of  the  numerical  one.  Figure  1  shows  the  wavenumbers  associated  with  the  Taylor  and  Kolmogorov 
scales  both  in  product  and  fuel  in  a  stoichiometric  H2-air  mixture.  These  values  were  determined  using  the  temperature 
dependence  of  the  viscosity  coefficient  similar  to  that  of  other  diffusion  coefficients,  namely  v  =  vo7’"/p  where  n  -  0.7 
and  Vo  =  2.9  x  10“^  g/s-cm-K°^  corresponding  io  Pr  -  Sc  -  1  [17].  It  can  be  seen  that  resolving  the  Kolmogorov 
scale  in  the  cold  fuel  would  require  substantially  higher  grid  resolutions.  Achieving  such  resolution,  while  maintaining 
reasonably  large  separation  between  the  scales  L  and  6l,  would  be  difficult  to  impossible  with  current  computational 
infrastructure.  For  higher  turbulent  intensities,  rjf  would  be  even  smaller.  In  addition,  physical  viscosity  results  in 
a  much  shallower  and  wider  dissipation  range.  Given  the  limited  range  of  spatial  scales  accessible  in  a  numerical 
simulation,  this  would  substantially  restrict  the  possible  extent  of  the  inertial  range  in  the  flow. 

Figure  1  shows  that  in  cases  SI  and  S3,  the  numerical  Taylor  scales,  i.e.,  the  largest  scales  in  the  domain  affected 
by  dissipation,  are  similar  to  physical  Taylor  scales  in  product  and  fuel,  respectively.  Therefore,  performed  simulations 
do  capture  both  the  smallest  and  the  largest  extents  of  the  inertial  range  present  in  the  physical  stoichiometric  H2-air 
mixture. 

These  considerations  show  that  fully  resolved  DNS  models  with  physical  viscosity  are  impractical  for  our  study. 
More  important,  however,  is  that  it  is  possible  to  perform  a  detailed  analysis  of  the  effects  of  small  scales  on  the 
turbulent  flame  by  varying  numerical  resolution  without  resorting  to  fully  resolved  DNS  simulations.  In  particular, 
if  a  certain  property  of  the  turbulent  flame  does  not  change  or  appears  to  be  converged  with  increasing  resolution,  it 
indicates  that  its  evolution  is  determined  only  by  large  scales.  On  average,  the  flow  on  large  scales  is  independent 
of  small  scales.  Therefore,  if  large-scale  motions  are  properly  captured  in  a  simulation,  as  is  the  case  in  the  models, 
presented  here,  then  the  evolution  of  such  a  property  cannot  be  affected  by  the  even  smaller  scales  that  would  exist  in 
a  real  system  due  to  a  much  smaller  physical  viscosity. 

Methods,  which  rely  on  numerical  viscosity  to  provide  kinetic-energy  dissipation  in  the  grid  while  explicitly 
capturing  the  inertial  range  of  the  energy  cascade,  are  often  referred  to  as  implicit  large-eddy  simulations  (ILES). 
Feasibility  of  this  approach  was  first  suggested  by  Boris  [28]  and  it  has  since  been  extensively  used  in  turbulence 
modeling.  In  particular,  a  detailed  analysis  of  the  effects  of  numerical  viscosity  in  the  context  of  the  PPM  method, 
similar  to  the  one  used  in  this  work,  was  performed  by  Porter  et  al.  [29].  An  example  of  similar  analysis  for  a  different 
higher-order  algorithm,  namely  flux-corrected  transport,  can  be  found  in  [30].  Comprehensive  overview  of  this  class 
of  methods  is  given  in  [31]. 

Aside  from  controlling  the  amount  of  energy  on  small  scales,  grid  resolution  also  plays  its  traditional  role.  It  must 
be  high  enough  to  provide  sufficient  accuracy  of  the  overall  solution  on  all  scales.  As  was  demonstrated  in  [19],  in 
Athena-RFX,  a  resolution  of  even  8  cells  per  6l  is  sufficient  to  reproduce  very  accurately  properties  of  the  planar 
laminar  flame.  A  priori,  it  is  not  clear,  however,  that  this  resolution  is  enough  to  capture  the  complex  dynamics  of 
the  turbulent  flame  which  contains  highly  distorted  flamelets.  Strong  turbulence  can  fold  individual  flame  sheets  with 
the  curvature  radius  comparable  to  6l  or  bring  two  flame  sheets  together.  This  would  form  cusps,  whose  velocity  of 
propagation  can  be  substantially  larger  than  that  of  a  planar  laminar  flame.  Properly  capturing  the  rate  of  burning 
in  such  cusps  depends  sensitively  on  the  code’s  ability  to  accurately  model  thermal  flux  focusing  in  regions  of  high 
flame  curvature.  Even  the  high-order  dimensionally  unsplit  algorithms  introduce  some  degree  of  anisotropy  to  the 
thermal  flux,  which  tends  to  align  itself  with  the  grid.  At  low  resolution,  this  can  substantially  affect  thermal  transport 
in  cusps.  Increasing  resolution  in  simulations  SI  -  S3  allows  us  to  carry  out  the  convergence  study  to  address  these 
issues. 

All  simulations  were  performed  in  the  domain  with  the  high  length-to-width  ratio  of  16: 1 .  The  longest  dimension 
of  the  domain  is  assumed  to  be  along  the  z-axis.  Such  long  domains  are  computationally  expensive,  however  they  are 
needed  to  be  able  to  follow  the  flame  evolution  for  extended  periods  of  time,  namely  for  IbTed-  We  learned  empirically 
that  this  extended  domain  size  is  necessary  to  accommodate  flame-brush  motion  with  respect  to  the  grid  during  that 
time  and  to  ensure  that  the  flame  brush  remains  sufficiently  far  from  the  open  outflow  boundaries  in  order  to  minimize 
their  effect  on  the  system. 

The  domain  in  all  calculations  was  initialized  with  uniform  density  po  and  temperature  To  (see  table  1).  In 
simulations  SI  and  S2,  initial  fluid  velocities  were  set  to  0.  In  contrast,  in  S3  the  velocity  field  was  initialized  with  the 
ideal  energy  spectrum  &{k)  oc  ^  extending  from  the  energy  injection  scale  L  to  the  numerical  Kolmogorov  scale 
7}  -  2Ax.  This  initial  spectrum  was  normalized  to  ensure  that  at  f  =  0  the  total  kinetic  energy  in  the  domain  was  equal 
to  its  predicted  steady-state  value  as  described  in  [19].  The  main  advantage  of  this  method  of  the  velocity  initialization 
is  that  it  allows  the  flow  to  reach  its  equilibrium  state  much  faster.  In  particular,  the  equilibration  time  is  2Teq  in  S3, 
whereas  it  is  3Teq  in  SI  and  S2  in  which  initial  velocities  are  zero  [19].  That  is  due  to  a  much  shorter  time  needed  to 
populate  all  spectral  modes  and  for  the  steady  energy  cascade  to  develop.  Moreover,  the  rate  of  energy  dissipation  is 
close  to  its  equilibrium  value  practically  from  f  =  0,  thus  preventing  the  build-up  of  excess  kinetic  energy  and,  thereby, 
also  shortening  the  equilibration  phase.  A  detailed  discussion  of  this  method  of  turbulent  flow  field  initialization  can 
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be  found  in  [19]. 

It  is  shown  in  [19]  that  in  nonreactive  turbulence,  the  long-term  system  evolution  is  not  sensitive  to  the  form  of 
the  velocity  field  in  the  domain  at  f  =  0,  provided  that  sufficient  time  was  given  for  the  system  to  reach  its  equilibrium 
state.  Further  results  will  demonstrate  that  this  is  also  the  case  in  reactive  turbulence. 

Steady  driving  with  constant  rate  e  of  energy  injection  into  the  domain  was  applied  to  the  system  for  the  total 
duration  of  all  simulations.  Energy  injection  scale  is  always  assumed  to  be  the  domain  width  L.  The  value  of  s 
was  chosen  to  produce  turbulence  field  of  sufficiently  high  intensity  which,  however,  was  low  enough  to  minimize 
the  probability  of  creating  the  weak  transonic  shocklets  that  arise  from  intermittency  in  the  turbulent  flow.  In  fast 
turbulent  flows,  which  can  still  be  nominally  characterized  as  subsonic  based,  for  instance,  on  the  value  of  Urms,  such 
features  can  represent  a  substantial  part  of  the  flow  field  [19]. 

Table  2  lists  a  number  of  velocity  characteristics  along  with  the  integral  scale  of  the  steady-state  turbulent  flow 
prior  to  the  moment  of  ignition.  In  the  simulations,  we  do  not  follow  the  evolution  of  the  nonreactive  turbulence  in 
the  steady  state  long  enough  to  extract  time-averaged  values  of  those  parameters.  Therefore,  table  2  lists  theoretically 
predicted  equilibrium  values  based  on  the  expressions  given  in  [19].  Analysis  performed  in  [19]  shows  that  such 
estimates  generally  have  accuracy  of  a  few  percent.  All  derived  quantities  listed  in  table  2,  such  as  Lc,  Da,  Maf, 
and  Map,  are  also  based  on  those  values. 

Given  the  large-scale  velocities  and  system  size  (or  integral  scale),  performed  simulations  can  be  characterized  as 
being  in  the  well-stirred  reactor  regime,  based  on  the  traditional  classification  of  [1,  2,  3].  In  particular,  the  Damkohler 
number.  Da,  defined  as  the  ratio  of  the  turbulent  integral  time  scale,  IjUs,  to  the  characteristic  chemical  time  scale, 
^lIS L,  is  0.1.  In  this  regime,  turbulent  transport  is  expected  to  overwhelm  molecular  diffusive  processes,  thereby 
significantly  affecting  the  internal  flame  structure  and  resulting  in  the  formation  of  broad  distributed  reaction  zones. 
The  Gibson  scale  Lc,  defined  as  the  scale  on  which  characteristic  turbulent  velocity  equals  S p,  is  almost  four  orders 
of  magnitude  smaller  than  dp.  This  suggests  that  turbulent  cascade  should  penetrate  deep  inside  individual  flamelets 
disrupting  them.  Performed  simulations  allow  us  to  investigate  the  validity  of  such  classification  of  this  combustion 
regime. 

In  our  model  we  do  not  include  any  artificial  cooling  mechanisms  to  compensate  for  the  turbulent  heating  of  the 
flow.  Based  on  the  analysis  given  in  [19],  the  amount  of  energy  injected  into  the  domain  within  one  large-scale  eddy 
turnover  time  Ted  is  ~  80%  of  the  total  steady-state  kinetic  energy  in  the  nonreactive  turbulent  flow.  The  typical  time 
step  is  ~  10“^  -  10^^  s,  or  ~  (see  table  2).  Therefore,  the  amount  of  energy  injected  into  the  domain  on  every 

time  step  is  <  10^“*  of  the  equilibrium  kinetic  energy. 

The  corresponding  relative  increase  in  the  internal  energy  and  temperature  due  to  the  turbulent  heating  in  the 
nonreactive  flow  within  Ted  is  [19] 


Et  -  E,fi 

Et,o 


T-Tq 

To 


^PO^ed 

E,fl 


-  Dy(y-  l)Map, 


(9) 


where  constant  D  -  0.5  and  Epo  is  the  thermal  energy  in  the  domain  at  f  =  0.  Using  values  of  y  from  table  1  and  Map 
from  table  2,  this  increase  is  ss  0.6%  of  their  initial  values  or,  equivalently,  ss  1.8  K  over  Ted-  Such  a  small  amount  of 
heating  cannot  result  in  any  significant  fuel  preconditioning. 

In  simulations  SI  and  S2,  the  flow  field  was  allowed  to  evolve  for  the  time  =  Med,  and  in  S3  for  the  time 
tign  -  2Ted  to  develop  the  steady-state  turbulent  flow  field.  At  the  time  f,g„,  a  planar  laminar  flame  with  its  front 
parallel  to  the  x  -  y  plane  was  imposed  in  the  domain.  The  initial  flame  position  zt,o,  given  in  table  2,  was  assumed 
to  be  the  location  of  the  Y  -  0.5  point  in  the  exact  laminar  flame  solution  for  the  reaction  model  parameters  and  fuel 
conditions  used.  Values  of  p,  P,  and  Y  in  each  cell  were  reset  to  those  obtained  from  the  laminar  flame  structure 
based  on  the  cell-center  z-coordinate.  Product  is  located  toward  the  left  z-boundary  and  fuel  is  located  toward  the 
right  z-boundary.  The  velocity  field  was  not  modified,  thereby  preserving  the  structure  that  developed  during  the 
equilibration  stage.  Hereafter,  for  simplicity  we  refer  to  the  moment  of  ignition  as  f  =  0. 

Prior  to  tign,  all  domain  boundaries  are  periodic.  At  tign,  boundary  conditions  (BCs)  along  the  left  and  right 
z-boundaries  are  switched  to  zero-order  extrapolation.  As  the  flow  evolves,  it  develops  complex  patterns  near  the 
boundaries.  Local  turbulent  motions  may  have  velocity  vectors  directed  both  in  and  out  of  the  domain  at  various  points 
along  the  boundary.  Moreover,  this  is  superimposed  on  the  large-scale  net  flow  associated  with  the  expansion  of  burnt 
material.  Consequently,  BCs  must  provide  the  possibility  for  the  fluid  to  move  both  in  and  out  of  the  domain  while 
adjusting  to  the  flow.  We  find  that  rather  simple  zero-order  extrapolation  BCs  successfully  manage  to  accomplish  this 
goal  while  preventing  any  unphysical  pressure  build-up  in  the  domain. 

This  type  of  BCs  is  known  to  cause  reflections  of  acoustic  perturbations.  High-speed  turbulent  flow  in  the  domain, 
however,  itself  generates  substantial  acoustic  noise  and  even  weak  transonic  shocklets.  Such  features  are  substan¬ 
tially  stronger  than  reflections  from  the  boundaries,  so  that  it  becomes  nearly  impossible  to  distinguish  such  reflected 
perturbations  from  the  ones  generated  by  the  flow  field  itself. 
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Finally,  fluid  entering  the  domain  with  the  large-scale  inflow  does  not  have  the  proper  K41-type  spectrum.  Pro¬ 
vided  that  the  turbulent  flame  brush  is  located  sufficiently  far  from  the  boundaries,  the  typical  time  for  the  fluid  to 
travel  from  the  boundary  to  the  flame  is  longer  than  a  few  Ted-  Therefore,  this  time  is  sufficient  for  the  incoming  flow 
to  be  subjected  to  the  external  driving  in  the  domain  and  to  develop  the  equilibrium  spectrum. 


3.  Results 


3.1.  Overall  evolution  and  global  properties  of  the  flame  brush 

Upon  ignition,  the  planar  laminar  flame  starts  being  wrinkled  by  turbulent  motions,  and  the  turbulent  flame  brush 
gradually  develops.  Typically,  within  ~  2Ted,  the  flame  brush  reaches  a  quasi-steady  state  in  which  its  width  6t  and 
speed  5  7  on  average  remain  constant. 

Figure  2  shows  the  three-dimensional  (3D)  structure  of  the  turbulent  flame  brush  based  on  the  distribution  of  the 
fuel  mass  fraction.  A  general  examination  of  the  flame  morphology  in  fig.  2  already  suggests  crucial  similarities 
as  well  as  differences  in  the  system  at  different  resolutions.  In  all  three  simulations,  the  flame  brush  represents  a 
highly  convolved  flame  with  a  thinner  reaction  zone  and  a  thicker  preheat  zone.  Thin  black  line  on  the  sides  of  the 
computational  domain  marks  the  boundary  between  the  preheat  and  reaction  zones,  and  the  thin  white  line  shows  the 
location  of  the  peak  reaction  rate.  Typically,  these  two  lines  closely  follow  each  other.  The  preheat  zone,  shown  in 
green  and  blue,  is  much  wider  than  the  reaction  zone,  and  generally  its  shape  does  not  follow  that  of  the  reaction  zone 
(see  also  fig.  4). 

The  flame  surface,  when  viewed  from  the  product  side,  appears  remarkably  independent  of  the  resolution.  In  all 
cases  it  is  smooth  on  smaller  scales  and  it  is  curved  on  larger  scales  comparable  to  the  turbulent  integral  scale.  In 
contrast,  the  flame  surface  on  the  fuel  side  is  wrinkled  on  progressively  smaller  scales  with  increasing  resolution.  This 
suggests  that  turbulent  motions  on  smaller  scales  affect  the  flame  structure  in  the  preheat  zone  and  their  effect  becomes 
less  pronounced  with  increasing  temperature  toward  the  reaction  zone. 

Figure  3  shows  the  evolution  of  6t  and  St,  normalized  by  6l  and  Sl-  Hereafter,  we  define  the  width  of  the 
turbulent  flame  brush  as 

dj  —  ^X.max  Zo.mm,  (10) 


where  zo.min  is  such  a  coordinate  that  Y  <  0.05  for  all  z  <  Zo.min  and  zi.max  is  such  that  Y  >  0.95  for  all  z  >  Zi.max- 
Therefore,  zo.min  marks  the  z-plane  to  the  left  of  which  is  pure  product,  and  zi.max  marks  the  z-plane  to  the  right  of 
which  is  pure  fuel.  This  is  illustrated  in  fig.  4. 

There  are  several  different  ways  to  define  the  turbulent  burning  velocity,  as  described  in  detail  in  [5].  The  flame 
brush  propagates  into  the  fuel  and  its  total  net  displacement  velocity  is  a  function  of  both  the  global  fuel-consumption 
speed  as  well  as  the  net  fuel  velocity  in  the  domain.  The  net  fuel  velocity  depends  sensitively  on  many  factors,  such  as 
the  details  of  the  flow  at  the  moment  of  ignition,  subsequent  evolution  of  the  stochastic  turbulent  field,  and  conditions 
at  the  domain  boundary.  Consequently,  the  speed  of  the  flame-brush  displacement  with  respect  to  the  grid  can  vary, 
and  it  cannot  be  predicted  a  priori. 

Since  we  are  considering  an  ideally  infinite,  unbounded  system,  there  are  no  boundaries  or  obstacles  with  respect 
to  which  the  velocity  of  the  flame  brush  could  be  defined.  At  the  same  time,  displacement  velocity  with  respect  to  the 
laboratory  reference  frame,  associated  with  the  computational  mesh,  does  not  have  physical  significance.  In  principle, 
we  could  associate  a  stationary  reference  frame  with  the  fuel  located  (infinitely)  far  from  the  flame  brush.  Such  fuel, 
while  being  in  the  state  of  turbulent  motion,  could  be  considered  stationary  on  average,  and  the  displacement  velocity 
of  the  flame  brush  with  respect  to  it  could  be  defined.  In  practice,  however,  the  volume  of  fuel  between  the  flame 
brush  and  the  domain  boundary  is  not  sufficiently  large  to  allow  us  to  perform  averaging  that  would  eliminate  spatial 
and  temporal  variability  in  the  fuel  velocity  and  would  provide  a  reliable  inertial  reference  frame. 

Therefore,  we  define  the  turbulent  flame  brush  velocity  5  7  as  the  global  fuel  consumption  speed  (cf.  eq.  (15)  in 


[5]) 


(11) 


Here  is  the  total  rate  of  fuel  consumption  inside  the  flame  brush,  i.e.,  the  total  mass  of  reactants  which  is  trans¬ 
formed  to  product  per  unit  time.  On  average  the  flame  brush  is  planar  in  the  domain  with  the  cross-sectional  area  of  L^. 
Then  in  steady  state  in  order  for  the  flame  brush  to  support  fuel  consumption  rate  hiR,  the  fuel  must  be  supplied  to  the 
flame  through  the  area  at  the  rate  poS  tL^-  Note,  that  even  though  the  fuel  density  varies  locally  by  a  small  amount 
due  to  the  compressibility  in  the  flow,  on  average  it  remains  equal  to  the  original  value  po.  Therefore,  in  a  perfectly 
steady  state,  S  7  represents  both  the  fuel  consumption  speed  as  well  the  velocity  of  the  flame  brush  displacement  with 
respect  to  the  stationary  fuel  far  from  the  flame  brush. 


Table  3  lists  time-averaged  values  of  St  and  S  t-  Time  averaging  is  performed  over  14  large-scale  eddy  turnover 
times.  The  start  of  the  averaging  interval  is  chosen  as  f  =  2Ted,  when  the  system  enters  the  quasi-steady  state.  While 
the  particular  moment  of  the  onset  of  this  state  is  not  precisely  defined,  we  find  that  the  flame-brush  parameters, 
discussed  in  this  paper,  reach  their  time-averaged  values  by  2Ted- 

Even  the  relatively  high-speed  turbulence,  considered  here,  only  increases  S  t  rather  modestly  with  the  value  of 
St  on  average  saturating  12m/s.  The  flame-brush  width  in  S3  is  somewhat  larger  than  the  domain  width, 

namely  St  -  1.8L,  and  it  is  7.7  times  larger  than  the  turbulent  integral  scale  1. 

Table  3  also  lists  the  order  of  self-convergence  of  each  parameter.  Since  the  computational  cell  size  decreases  by 
a  factor  of  2  progressively  between  each  simulation,  the  order  of  self-convergence  (9(^6)  of  a  parameter  (p  is  defined  as 


0(0)  =  log2 


I0r»«l  -  0r««3l 

J0ran2  -  0ran3l^ 


(12) 


The  table  shows  that  both  St  and  S  t  converge  quadratically,  as  would  be  expected  for  a  second-order  numerical 
method. 

In  a  dynamic  unsteady  flow,  fuel  consumption  inside  the  flame  brush  is  not  perfectly  balanced  by  the  influx  of 
fresh  fuel.  Inflow  of  fuel  can  dominate  as  the  flame  brush  increases  its  width  and  the  amount  of  fuel  inside.  This  is 
followed  by  a  period  when  fuel  consumption  prevails,  rapidly  burning  the  accumulated  reactants  and  causing  the  flame 
brush  to  shrink  in  size.  While  on  average  those  two  processes  do  remain  balanced  maintaining  constant  average  flame 
brush  width  and  speed,  as  was  just  discussed,  such  variability  is  an  important  part  of  the  overall  flame-brush  evolution. 
Figure  3  shows  that  in  the  lowest  resolution  case  SI,  bj  varies  by  more  than  a  factor  of  two,  and  S t  varies  by  more 
than  a  factor  of  three.  Furthermore,  the  overall  variability  decreases  with  increasing  resolution.  This  is  particularly 
apparent  for  the  turbulent  flame  speed.  Note  that  peak  values  of  5 7  in  S3  are  <6,  while  in  SI  they  reach  ss  14  in  one 
episode. 

Figure  2,  which  shows  the  flame  brush  structure  in  all  three  simulations  at  t  =  \2>Ted,  also  illustrates  such  variable 
nature  of  the  system  evolution.  Figure  3  shows  that  at  this  time,  SI  and  S3  undergo  two  extreme  episodes  in  their 
evolution.  The  flame  brush  in  SI  reaches  its  maximum  width,  which  is  ~  50%  larger  than  its  average  value,  while  S  7 
increases  by  more  than  a  factor  of  two.  S3,  on  the  other  hand,  is  in  the  most  quiescent  phase  with  St  shrinking  by 
one  third  and  S  7  being  only  twice  the  laminar  flame  speed.  Finally,  S2  at  this  time  is  in  its  average  state  with  both 
St  and  S  7  very  close  to  their  average  values.  Therefore,  aside  from  the  differences  in  the  preheat  zone  structure,  the 
three  states  of  the  turbulent  flame  shown  in  fig.  2  represent  the  main  stages  in  the  flame-brush  evolution  which  are 
characteristic  of  all  three  simulations.  Namely,  they  illustrate  recurring  transitions  between  the  periods  of  a  widened 
flame  brush,  containing  highly  convolved  flame  surface,  and  more  quiescent  stages,  when  the  thinner  flame  brush 
contains  a  smoother,  flatter  flame. 

Qualitatively,  the  nature  of  such  variability  can  be  explained  as  follows.  The  internal  structure  of  the  flame  brush  is 
determined  by  two  counteracting  processes.  On  one  hand,  flame  surface  area  is  created  by  turbulent  motions  bending 
and  stretching  the  flame  and,  thereby,  increasing  both  the  flame  brush  thickness  and  the  flame  surface  density  inside 
the  brush.  On  the  other  hand,  this  process  is  balanced  by  flame-surface  destruction,  which  is  an  inherently  nonlinear 
stabilization  process  preventing  unbounded  growth  of  the  flame  perturbations  [32,  33].  Flame-surface  destruction 
takes  place  in  regions  where  the  curvature  radius  of  the  flame  surface  becomes  comparable  to  the  flame  thickness 
(“cusps”)  or  where  the  flame  sheets  come  close  to  each  other.  Such  a  picture  of  the  balance  between  surface  creation 
and  destruction  has  been  reflected  in  a  number  of  analytical  and  numerical  models  of  turbulent  flames  (e.g.,  see 
[34,  35]  as  well  as  the  review  by  [5]  and  references  therein). 

Three  distinct  stages  of  the  destruction  process  can  be  seen  in  the  upper  and  middle  panels  of  fig.  2.  Consider 
region  A,  where  the  flame  surface  is  tightly  packed  by  turbulence.  This  results  in  two  flame  sheets  being  so  close  to 
each  other  that  their  preheat  zones  start  to  overlap.  Moreover,  each  flamelet  is  also  curved  on  scales  comparable  to  its 
thickness.  Thermal  flux  from  hot  products  into  the  colder  fuel  becomes  focused  in  such  cusps,  thus  accelerating  the 
heating  of  reactants  and  effectively  increasing  the  local  burning  velocity.  This  makes  the  two  flame  sheets  propagate 
faster  toward  each  other. 

Eventually,  this  may  lead  to  the  configuration  seen  in  region  B.  There  preheat  zones  of  each  flame  sheet  substan¬ 
tially  overlap  even  though  the  reaction  zones  are  still  separated.  At  this  point,  the  whole  region  acts  as  one  large  cusp 
with  the  heat  flux  coming  from  all  directions  and  thus  rapidly  heating  up  the  unreacted  material.  The  temperature 
quickly  reaches  the  ignition  point  and  all  of  this  region  becomes  a  reaction  zone. 

Examples  of  such  merged,  extended  reaction  zones  can  be  seen  in  regions  C.  Those  parts  of  the  flame  quickly 
bum  out,  which  reduces  the  curvature  and  the  overall  surface  area  of  the  flame.  Such  events  are  transient  phenomena, 
and  they  must  not  be  mistaken  for  stable  distributed  reaction  zones.  In  particular,  it  is  very  important  in  experimental 
settings  for  the  diagnostics  to  be  able  to  properly  characterize  such  broad  reaction  zones  as  transients. 
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In  the  fast  turbulence,  as  present  in  these  simulations  with  Da  1,  the  characteristic  turbulent  time  dilUg  at 
the  scale  of  the  laminar  flame  width  is  much  smaller  than  the  laminar  flame  self-crossing  time  SlIS l-  Therefore, 
turbulence  can  bring  two  flame  sheets  together  and  then  pull  them  apart  before  each  one  propagates  over  the  distance 
~  6l-  As  a  result,  flame  sheets  in  region  A  in  fig.  2  can  be  separated  by  turbulence  before  they  ever  merge.  On  the 
other  hand,  the  local  burning  velocity  in  cusps  can  be  very  large,  in  fact,  ideally  it  can  be  infinite.  Therefore,  once 
a  cusp  is  formed,  i.e.,  there  exists  a  region  of  significant  flux  focusing,  it  burns  out  on  time  scales  short  enough  for 
turbulence  essentially  to  have  no  effect  on  it. 

In  summary,  periodically  the  flame  brush  undergoes  a  state  in  its  evolution  when  turbulence  rapidly  increases  the 
flame  surface  which  becomes  highly  convolved  and  develops  cusps.  This  situation  is  represented  by  the  upper  panels 
in  fig.  2.  The  local  burning  velocity  in  cusps  increases  substantially  to  the  point  when  turbulent  motions  are  no  longer 
locally  dominant.  This  causes  cusps  to  bum  out  quickly  leading  to  an  increased  rate  of  global  fuel  consumption  and 
smoothing  of  the  overall  flame  surface.  This  stage  is  shown  in  the  middle  panels  in  fig.  2.  Eventually,  this  results 
in  a  flame  brush  which  is  substantially  thinner  and  which  contains  much  less  convolved  flame  surface,  as  seen  in  the 
bottom  panels  of  fig.  2.  Such  flatter,  slower  flame  with  fewer  cusps  now  again  becomes  susceptible  to  the  action  of 
turbulent  motions,  which  can  increase  its  surface  area  and,  hence,  the  cycle  repeats. 

Finally,  fig.  5  shows  temperature  structure  in  the  flame  bmsh  in  S2,  corresponding  to  the  middle  panels  of  fig.  2. 
As  expected,  the  temperature  distribution  closely  follows  that  of  the  fuel  mass  fraction. 

3.2.  Internal  structure  of  the  flame  brush 

Figure  2  suggests  that  even  in  the  presence  of  high-speed  turbulence,  the  flame  brush  consists  of  highly  convolved 
flamelets.  In  order  to  examine  this,  we  need  to  identify  such  flamelets  and  determine  their  internal  stmcture,  if  they 
are  indeed  present.  To  accomplish  this,  we  developed  the  following  method. 

The  range  of  values  of  Y  and  T  in  the  domain  is  discretized  into  n  equal  intervals.  Typically,  we  choose  Y  e 
[0.01, 0.99]  and  T  6  [350.0, 2135.0]  K.  The  upper  bound  of  the  temperature  interval  is  the  adiabatic  flame  temperature 
Tp.  The  lower  bound  is  somewhat  higher  than  To.  In  the  course  of  the  simulation,  fuel  temperature  away  from  the 
flame  slowly  rises  due  to  the  turbulent  energy  dissipation,  thereby  raising  the  minimum  temperature  value  in  the 
domain  (see  §  2.4).  The  lower  bound  of  the  temperature  interval  is  chosen  to  account  for  the  heating  of  the  fuel  and 
to  include  only  the  temperature  increase  associated  with  flame  propagation.  For  each  of  the  discrete  values  of  T,  and 
r„  with  i  taking  values  0  to  n,  we  construct  an  isosurface  Si  with  area  A,  using  the  “marching  cubes”  algorithm.  We 
also  determine  the  volume  V,  bounded  by  the  isosurfaces  Si  and  >S,+i.  Since  the  flame  brush  always  intersects  the 
computational  domain  boundaries,  this  volume  will  always  be  bounded  by  Si,  >S,+i,  and  the  boundaries.  Then  we 
construct  the  expression 

mn  =  VilAi,  (13) 

assuming  rjo  -  0.  The  rji+i  gives  the  average  distance  separating  all  points  in  the  domain  in  which  Y  —  T,+i  and 
T  -  r,+i  from  those  in  which  Y  -  Yi  and  T  -  Ti  respectively.  Finally,  Yirf)  and  Tirf)  give  the  average  flame  let 
structure  in  the  flame  brush.  In  Appendix  A  we  discuss  in  detail  the  properties  and  limitations  of  this  method,  as  well 
as  the  issues  associated  with  its  practical  implementation. 

For  finite  values  of  AT  s  T,+i  -  T,  or  AT  s  Ti+i  -  Ti,  the  local  separation  between  the  two  isosurfaces  Si  and 
>S,+i  will  vary  from  point  to  point.  Consider  fig.  2.  Thin  white  and  black  lines  mark  isosurfaces  T  =  0.2  and  T  =  0.6. 
Regions  can  be  identified  where  they  come  very  close  to  each  other  and  where  they  are  substantially  separated,  e.g.,  as 
in  regions  C  (see  also  fig.  4).  The  same  is  also  true  for  temperature,  as  evidenced  by  fig.  5.  We  show  in  Appendix  A, 
however,  that  as  AT  and  AT  become  smaller,  the  variation  in  the  local  separation  of  the  two  isosurfaces  will  also 
decrease. 

Therefore,  rji+i,  as  defined  in  eq.  (13),  can  be  considered  as  a  distance  between  consecutive  isosurfaces  averaged 
over  their  entire  area.  Regions  of  large  flame  stretch,  in  which  the  flamelet  becomes  very  thin,  would  tend  to  decrease 
the  value  of  rf.  On  the  other  hand,  cusps,  regions  of  significant  flame  folding,  and  broadened  reaction  zones  would 
tend  to  make  tj  larger.  Thus,  the  flame  structure  Y(ri)  and  T(ri),  produced  by  this  analysis  for  a  given  instantaneous 
state  of  the  flame  brush,  is  an  average  of  all  local  realizations  of  the  flame  structure  at  each  point  of  the  flame  surface. 

Subsequently,  Y(ti)  and  T{r])  can  be  time-averaged,  as  a  substitute  for  ensemble-averaging,  to  obtain  the  statis¬ 
tically  dominant  structure  of  the  flamelets  inside  the  flame  brush.  The  resulting  structure  then  contains  information 
about  the  statistical  weight  of  cusps  and  broadened  reaction  zones  in  the  ensemble  of  all  possible  flamelet  configura¬ 
tions  that  can  exist  for  given  conditions.  Similarity  of  this  structure  to  that  of  the  planar  laminar  flame  would  indicate 
that  cusps  and  broadened  reaction  zones  are  indeed  transient  phenomena,  and  they  do  not  dominate  the  internal  flame¬ 
brush  structure.  On  the  other  hand,  if  substantial  flame  broadening  by  turbulent  transport  indeed  takes  place,  as  has 
been  postulated  for  such  intense  turbulence,  then  broadened  reaction  zones  would  dominate  the  statistical  ensemble 
of  possible  flame  configurations,  and  the  resulting  averaged  structure  would  be  significantly  wider  than  the  laminar 
flame  profile. 


10 


We  applied  this  method  to  the  flow  in  simulations  SI  -  S3.  For  SI  and  S2,  n  -  6lI ^x,  as  given  in  table  2,  in 
accordance  with  the  criterion  in  eq.  (19).  For  S3,  however,  n  was  limited  to  a  lower  value,  namely  16,  due  to  the  high 
computational  cost  of  reconstructing  a  substantially  larger  number  of  isosurfaces'  in  the  case  of  n  =  32.  The  Y  and 
T  profiles  were  averaged  over  the  time  interval  (2  -  16)Ted-  In  particular,  averaging  was  performed  over  401  discrete 
time  states  in  SI,  101  in  S2,  and  419  in  S3.  The  resulting  profiles  are  shown  in  fig.  6. 

Values  of  the  x-coordinate  in  the  figure  are  chosen  such  that  y(x  =  0)  =  0.5  both  in  the  exact  laminar  flame  solution 
and  in  the  obtained  profiles.  In  the  exact  solution,  this  also  uniquely  determines  the  position  of  the  temperature  profile 
since  there  is  a  precise  relation  between  T  and  Y.  At  the  same  time,  the  method  we  are  using  to  reconstruct  the  flamelet 
structure  does  not  provide  any  means  to  relate  the  obtained  distributions  of  different  variables,  as  it  is  discussed  in 
detail  in  Appendix  A.  Values  of  rj  are  independent  in  each  profile.  It  is  not  a  common  spatial  coordinate,  and  it 
simply  indicates  distance  between  consecutive  isosurface  values  of  Y,  T,  etc.  Each  profile  is  a  statistically  averaged 
representation  of  all  existing  flamelet  configurations  in  the  flame  brush.  Thus,  any  relation  between  the  Y  and  T 
profiles  can  only  exist  in  the  statistical  sense,  namely  such  relation  would  only  indicate  the  most  likely  values  of  T  for 
given  values  of  Y  and  vice  versa.  Therefore,  here  for  simplicity  we  offset  temperature  profiles  so  that  the  coordinate 
of  the  temperature  value,  which  corresponds  to  the  reactedness  c  -  (T  -  To)l(Tp  -  Tq)  -  0.5,  would  be  x  =  0,  as  in 
the  exact  laminar  flame  solution. 

The  principal  conclusion  that  emerges  from  fig.  6  is  that,  on  average,  the  turbulent  flame  brush  interacting  with 
high-speed  turbulence  represents  highly  convolved  flamelets  that  have  internal  structure  of  a  laminar  flame  with  a 
somewhat  broadened  preheat  zone.  They  Y  and  T  profiles  are  very  close  to  the  exact  laminar  solution  in  the  reaction 
zone.  They  begin  to  diverge  from  the  laminar  profiles  only  at  the  outer  boundary  of  the  reaction  zone,  around  x/d/,  ss 
0.1  -  0.2.  In  fact,  the  Y  profile  for  S3  is  virtually  indistinguishable  from  the  laminar  solution  inside  the  reaction  zone. 
The  remarkable  similarity  of  the  flamelet  structure  inside  the  reaction  zone  to  that  of  the  planar  laminar  flame  is  a 
partial  justification  for  our  choice  of  the  x-coordinate  origin  for  the  temperature  profile. 

The  deviation  of  the  Y  and  T  profiles  from  the  laminar  solution  increases  with  distance  from  the  reaction  zone. 
Even  for  such  strong  turbulence,  however,  this  difference  is  fairly  small,  with  the  total  preheat  zone  width  being 
increased  by  less  than  a  factor  of  two.  Temperature  profiles  also  show  some  deviation  from  the  laminar  solution  near 
the  high-T  end.  This  is  mostly  due  to  errors  introduced  by  temperature  variations,  which  occur  away  from  the  flame 
brush  in  the  product  (see  Appendix  A  for  further  discussion  of  the  effects  of  such  variations  around  the  upper  bound 
of  the  discretized  temperature  range). 

Einally,  note  also  that  the  T  profiles  have  a  somewhat  higher  minimum  value  than  the  laminar  solution.  This  is  the 
consequence  of  fuel  heating  by  turbulence,  as  was  discussed  in  §  2.4. 

The  profiles,  obtained  for  all  three  simulations,  are  very  close  to  each  other  for  all  values  of  both  Y  and  T .  Given 
the  statistical  nature  of  the  distributions,  this  serves  as  evidence  that  the  simulations  can  be  considered  converged. 
More  importantly,  the  agreement  between  profiles  in  the  preheat  zone  shows  that  thermal  conduction  is  enhanced 
predominantly  by  turbulent  motions  on  scales  ~  6l-  The  amount  of  energy  contained  on  this  scale  is  the  same  in  all 
simulations  (see  fig.  1).  On  the  other  hand,  the  energy  contained  on  smaller  scales,  and,  thus,  the  velocity,  increases 
substantially  from  case  SI  to  S3.  If  those  scales  were  to  contribute  noticeably  to  the  overall  thermal  transport,  such 
increase  in  energy  content  on  those  scales  would  be  reflected  in  a  broader  preheat  zone.  There  is,  however,  no  evidence 
of  this. 

Distributions  of  both  Y  and  T  show  fairly  low  variability  throughout  the  time  interval  used  for  averaging.  In  the 
case  of  S3,  profiles  for  all  individual  time  states  are  contained  within  the  shaded  cyan  areas.  This  supports  our  choice 
of  this  interval  as  a  period  of  quasi-steady  state  in  the  system  evolution.  There  is,  however,  more  variability  in  the 
preheat  zone  than  in  the  reaction  zone.  Moreover,  this  is  more  prominent  for  fuel  mass  fraction  than  temperature. 

Einally,  it  is  important  to  emphasize,  that  in  the  preheat  zone  the  laminar  solutions  lie  outside  the  range  of  variabil¬ 
ity  of  the  obtained  profiles.  This  confirms  that  the  observed  broadening  is  statistically  significant,  and  such  flamelet 
structure  is  prevalent  throughout  the  evolution  of  the  system.  On  the  other  hand,  the  fact  that  the  laminar  solution  is  so 
close  to  the  boundary  of  the  variability  range  shows  how  weak  this  effect  is,  even  in  the  case  of  such  strong  turbulence 
as  is  present  in  these  simulations. 

4.  Discussion  and  Conclusions 

We  have  explored  the  dynamics  and  properties  of  a  turbulent  flame  formed  in  the  presence  of  subsonic,  high-speed, 
homogeneous,  isotropic  K41-type  turbulence  in  an  unconfined  system,  i.e.,  in  the  absence  of  walls  and  boundaries. 


’  Due  to  the  adaptive  nature  of  the  algorithm,  which  implements  the  reconstmction  method  as  described  in  Appendix  A,  the  actual  number  of 
isosurfaces  and  isovolumes  that  need  to  be  reconstructed  for  each  variable  can  be  significantly  larger  than  n.  In  particular,  for  n  =  32  on  average 
more  than  100  isosurfaces  need  to  be  reconstructed  for  each  time  state  with  their  total  number  exceeding  50,000. 
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A  Damkohler  number  for  the  turbulent  flow  is  Da  =  0.1  while  the  Gibson  scale  is  Lq  -  2.96  x  10  Therefore, 
the  system  considered  represents  turbulent  combustion  in  the  well-stirred  reactor  regime,  according  to  the  traditional 
classification  [1]. 

Numerical  modeling  was  performed  using  the  massively  parallel,  fully  compressible,  higher-order,  dimensionally 
unsplit,  reactive-flow  code  Athena-RFX  [19].  The  highest  resolution  calculation  S3  required  approximately  100,000 
CPU  hours  on  1000  processor  cores  of  the  Ranger  platform  at  TACC,  and  a  similar  amount  of  computational  time  was 
necessary  for  data  processing  and  reconstruction  of  the  flamelet  structure. 

4.7.  General  properties  of  the  turbulent  flame 

Our  results  show  that  turbulence-flame  interaction  leads  to  the  formation  of  a  steadily  propagating  turbulent  flame 
brush.  We  followed  the  system  evolution  in  this  quasi-steady  state  for  14  large-scale  eddy  turnover  times,  which  was 
sufficient  to  provide  long-term  statistics.  On  average,  both  the  flame  width  and  speed  remain  constant,  even  though  at 
each  moment  in  time  they  can  vary  substantially.  We  observed  an  increase  in  the  time-averaged  turbulent  flame  speed, 
St,  compared  to  the  laminar  flame  speed,  5/,,  by  a  factor  of  4.  The  time-averaged  turbulent  flame  width,  6t,  was 
almost  twice  the  domain  width  and  the  energy  injection  scale,  L.  It  was  also  almost  8  times  larger  than  the  turbulent 
integral  scale,  1. 

An  important  question  concerns  the  parameters  which  predominantly  determine  Sj-  In  our  simulations,  the  system 
size,  which  is  set  by  the  domain  width,  is  equal  to  L  and,  thereby,  it  is  close  to  1.  It  is,  therefore,  impossible  to  separate 
the  individual  role  played  by  the  system  size,  L,  and  I  in  defining  the  equilibrium  value  of  St-  In  particular,  given 
fixed  values  of  L  and  Z,  it  is  not  known  whether  a  substantial  increase  in  system  size  would  result  in  a  change  in  St- 
Moreover,  it  is  not  clear  whether,  in  such  a  large  system,  there  would  still  exist  the  observed  steady  flame  propagation 
with  a  fixed  width  which  is  much  smaller  than  the  system  size.  Finally,  an  important  question  concerns  the  relation 
between  L  and  Sj-  Does  an  increase  in  L  lead  to  a  corresponding  increase  in  6t,  or  is  there  a  limiting  value  of  djl 

The  study  by  Aspden  et  al.  [15],  which  used  faster  turbulence  relative  to  S l  than  this  work,  also  observed  a 
turbulent  flame  propagating  with  a  steady  width  and  speed.  In  that  work  energy  as  well  was  injected  at  the  scale  close 
to  the  domain  width.  In  their  fastest  calculation,  however,  the  domain  was  almost  8  times  larger  than  in  our  case 
relative  to  the  laminar  flame  width,  namely  ss  65(5l  vs.  8(5/,.  Yet  dj  was  ss  3  times  the  domain  width,  which  is 
only  50%  larger  than  in  our  simulations.  This  suggests  that  the  width  of  the  flame  brush  increases  in  a  way  similar  to 
the  increase  in  the  driving  scale  and  system  size.  Their  work,  however,  as  well  does  not  address  the  question  of  the 
individual  role  of  these  two  scales. 

In  systems  which  are  much  larger  than  the  driving  scale,  most  likely  both  L  and  I,  rather  than  the  domain  size, 
play  the  primary  role  in  setting  the  equilibrium  width  of  the  flame  brush.  Indeed,  the  turbulent  integral  scale  along 
with  the  driving  scale  are  the  characteristic  scales  of  coherent  turbulent  motions.  Moreover,  on  scales  T  »  L  the 
turbulent  energy  spectrum  has  the  dependence  [36].  Therefore,  the  energy  contained  in  those  large  scales  decays 
rapidly,  which  minimizes  their  potential  role  in  defining  the  turbulent  flame  width.  At  the  same  time,  in  situations 
when  domain  width  is  close  to  L,  domain  boundaries  imprint  artificial  periodicity  on  the  flow.  This  can  constrain  the 
growth  of  6t,  thereby,  setting  its  equilibrium  value. 

The  main  reason  for  the  difficulty  with  addressing  these  issues  is  the  limited  range  of  spatial  scales  presently 
accessible  even  in  the  most  heroic  numerical  calculations.  The  need  to  provide  sufficiently  high  resolution  at  the 
scale  of  the  laminar  flame,  and  still  maintain  some  separation  between  L  and  6l,  does  not  allow  the  system  size  to  be 
more  than  an  order  of  magnitude  larger  than  L.  This  is  insufficient  to  provide  a  conclusive  answer  to  the  questions 
mentioned  above. 

4.2.  Internal  structure  of  the  flame  brush  and  the  role  of  small-scale  turbulence 

Our  results  show  that  in  the  presence  of  high-speed  turbulence,  the  flame  brush  is  comprised  of  highly  convolved 
flamelets.  The  time-averaged  internal  structure  of  those  flamelets  is  very  close  to  that  of  the  planar  laminar  flame,  and, 
in  fact,  they  are  virtually  identical  in  the  reaction  zone.  The  preheat  zone,  however,  does  show  evidence  of  broadening, 
although  this  effect,  while  statistically  significant,  is  fairly  small  with  the  width  of  the  preheat  zone  increasing  by  less 
than  a  factor  of  two. 

The  internal  structure  of  flamelets  inside  the  flame  brush  was  reconstructed  directly  by  a  method  based  on  isosur¬ 
faces  of  the  3D  distribution  of  a  scalar  quantity,  such  as  fuel  mass  fraction,  Y,  or  temperature,  T .  The  isosurfaces  are 
obtained  at  evenly  spaced  values  of  Y  and  T  and  the  overall  method  does  not  make  any  assumptions  concerning  the 
possible  underlying  flame  structure  or  even  the  presence  of  the  flame  itself. 

Our  results  demonstrate  that  the  flamelet  reaction  zone  is  not  broadened  in  this  regime,  which  normally  would  be 
classified  as  a  well-stirred  reactor.  Turbulent  velocities  on  the  scale  of  the  laminar  flame  width  are  15  times  larger  than 
the  laminar  flame  speed.  In  principle,  this  would  suggest  that  the  turbulent  transport  should  overwhelm  the  molecular 
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thermal  conduction  and  species  diffusion.  It  appears,  however,  that  disrupting  the  flame  would  require  substantially 
higher  turbulent  intensities.  This  result  is  in  general  agreement  with  the  behavior  observed  by  [10,  15]. 

The  obtained  flamelet  structure  represents  the  statistically  dominant  state  of  the  flamelets  inside  the  flame  brush. 
At  the  same  time,  substantially  broadened  reaction  zones  can  arise  in  the  flow  as  a  result  of  collisions  and  mergers  of 
flame  sheets  or  in  the  regions  of  high  flame  curvature  (cusps).  Such  broadened  reaction  zones,  however,  are  transient 
in  nature  being  the  result  of  the  rapidly  changing  3D  configuration  of  the  flame  surface.  They  should  not  be  considered 
as  evidence  of  a  distributed  mode  of  burning  since  they  are  in  no  way  supported  by  the  turbulent  transport. 

The  progressive  increase  in  resolution  between  simulations  SI  and  S3  results  in  a  substantial  increase  in  the 
energy  of  turbulent  motions  on  scales  A  <  Sl  and  causes  the  energy  cascade  in  nonreactive  turbulence  to  extend 
to  much  smaller  scales.  This,  however,  is  insufficient  to  disrupt  the  internal  structure  of  the  reaction  zone  which 
is  demonstrated  by  the  fact  that  the  Y  and  T  profiles,  describing  the  internal  structure  of  the  flamelets  in  all  three 
simulations,  are  very  close  to  each  other  for  all  values  of  Y  and  T .  Furthermore,  since  the  structure  of  the  broadened 
preheat  zone  does  not  change  between  SI  and  S3,  this  suggests  that  such  broadening  is  determined  by  scales  d  >  <5/, 
rather  than  small-scale  motions.  The  effect  of  small-scale  turbulence  appears  only  in  the  progressively  finer  wrinkling 
of  the  flame  surface  on  the  fuel  side  in  the  coldest  part  of  the  preheat  zone. 

These  results  lead  to  the  following  important  conclusion.  The  turbulent  energy  cascade  fails  to  penetrate  the 
internal  structure  of  the  flame  on  scales  much  smaller  than  6l  even  in  the  presence  of  high-speed  turbulence  considered 
in  this  work.  This  shows  that  the  response  of  the  flame  to  the  action  of  turbulence  is  qualitatively  different  from  that  of 
a  passively  advected  scalar,  which  can  generally  be  a  source  of  only  limited  insight  into  the  behavior  of  the  turbulent 
flames. 

What  mechanisms  are  responsible  for  suppressing  the  effects  of  small-scale  motions  in  the  reactive  turbulence? 
Dramatic  difference  in  the  degree  of  wrinkling  of  the  flame  surface  on  the  fuel  and  product  sides,  which  was  observed 
in  the  highest  resolution  case  S3,  suggests  that  this  suppression  of  small-scale  turbulence  occurs  as  the  flow  passes 
through  the  individual  flamelets  rather  than  the  whole  flame  brush.  Two  possible  processes,  which  agree  with  this 
observation,  were  previously  suggested  (see  discussion  in  [5]  and  references  therein).  On  one  hand,  as  the  fluid  passes 
through  the  flame,  it  heats  up  and  expands.  Since  the  circulation  of  individual  eddies  must  be  conserved  in  the  process 
of  this  expansion,  the  rotational  velocity  of  eddies  will  significantly  decrease.  Moreover,  this  will  also  shift  energy 
from  smaller  to  larger  scales.  On  the  other  hand,  such  expansion  causes  rapid  acceleration  of  the  fluid  in  the  reference 
frame  of  the  flame  brush.  As  a  result,  the  residence  time  of  individual  eddies  substantially  decreases  as  the  fluid 
quickly  leaves  the  flame  brush. 

It  is  not  clear,  however,  whether  these  two  mechanisms  are  sufficient  to  completely  account  for  the  observed 
behavior.  It  is  also  not  known  whether  their  effect  on  small  and  large  scales  is  similar  or  not.  Answering  these 
questions  is  crucial  for  our  understanding  of  the  turbulence-flame  interaction. 

Ultimately,  our  understanding  of  the  turbulent  combustion  process  is  measured,  to  a  large  extent,  by  our  ability 
to  predict  the  turbulent  burning  velocity.  According  to  the  Damkohler  concept  [4]  discussed  in  §  1,  two  processes 
can  determine  the  turbulent  flame  speed;  increase  in  the  flame  surface  area  by  turbulent  stretching  and  folding  of  the 
flame,  and  increase  in  the  local  flame  speed  by  turbulent  diffusion. 

Similarity  of  the  internal  flame  structure  to  that  of  the  planar  laminar  flame  shows  that  local  flame  speed  cannot 
be  substantially  increased  by  turbulence.  The  flame  structure  is  determined  by  the  balance  of  the  reactions  and  the 
diffusive  processes.  In  our  model,  the  reaction  rate  is  determined  purely  by  the  local  thermodynamic  state  of  the  fluid 
and,  thereby,  it  is  not  directly  affected  by  turbulence  which  can  only  enhance  the  diffusive  transport.  For  a  given 
one-step  Arrhenius  reaction  model,  the  coefficients  of  thermal  conduction  and  species  diffusion  uniquely  determine 
both  the  flame  structure  and  speed.  If  turbulence  does  indeed  increase  the  effective  diffusive  coefficients,  this  would 
increase  the  local  flame  speed,  but  this  would  also  change  the  flame  structure.  We  do  not  find  any  evidence  of  this. 

This  raises  the  question  whether  the  observed  turbulent  flame  speed,  shown  in  fig.  3b,  can  be  fully  accounted  for 
only  by  the  increase  in  the  flame  surface  area  or  whether  other  processes  contribute  to  it.  We  will  discuss  this  issue  in 
detail  in  a  separate  paper  [37]. 

Finally,  it  is  important  to  emphasize  that  averaged  distributions  of  various  quantities,  such  as  Y,  T,  or  the  reaction 
rate,  R,  through  the  turbulent  flame  do  not  provide  information  about  the  internal  structure  of  the  flame  brush,  and 
they  cannot  be  viewed  as  evidence  of  flame  broadening.  Consider  fig.  7  which  shows  the  x-y-averaged  profiles  of  Y, 
T,  and  R  along  the  z-axis  in  simulation  S3.  All  profiles,  including  that  of  the  reaction  rate,  are  smooth  and  the  overall 
structure  is  similar  to  the  one  shown  in  fig.  6,  with  the  key  exception  that  the  width  of  the  profiles  is  substantially 
larger  by  a  factor  ss  4.  In  particular,  the  broad  smooth  profile  of  the  reaction  rate  can  lead  to  the  erroneous  conclusion 
that  this  system  represents  a  distributed  mode  of  burning  with  a  substantially  broadened  reaction  zone.  If  the  flame 
brush  volume  is  large  enough,  then  the  average  profiles  of  Y,  T,  etc.,  will  always  be  smooth  and  wide  since  they 
are  simply  the  result  of  averaging  many  individual  flamelet  structures.  Moreover,  the  turbulent  flame  is  the  region 
in  which  Y  and  T  transition  from  their  values  in  the  fuel  to  those  in  the  product.  Therefore,  those  profiles  will 
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smoothly  and  monotonically  change  from  their  fuel  to  product  values  and  the  averaged  distribution  of  R  will  simply 
vary  accordingly. 

4.3.  Convergence  and  resolution 

What  do  our  results  say  about  the  minimum  resolution  needed  to  provide  a  converged  solution  in  simulations  of 
high-speed  turbulent  combustion?  Based  on  table  3,  time-averaged  values  of  6tI6l  and  StISl  converge  quadrati- 
cally.  The  values  of  6tI6l  quantities  vary  only  by  ss  3%  in  the  two  highest-resolution  simulations  S2  and  S3,  while 
StISl  varies  by  10%.  Profiles  of  Y  and  T,  which  represent  the  internal  flamelet  structure,  also  demonstrate  excellent 
convergence. 

The  only  exception  to  this  converged  behavior  is  the  degree  of  flame  surface  wrinkling  on  the  fuel  side.  In  fact, 
convergence  of  this  characteristic  of  the  turbulent  flame  was  not  expected  at  the  resolutions  used  and,  indeed,  it  was 
not  observed.  It  was  discussed  in  §  2.4  that  simulations  with  different  computational  cell  sizes  do  not  represent  the 
same  physical  system.  The  spectral  energy  distribution  in  the  domain  changes  with  increasing  resolution  because  the 
amount  of  energy  contained  on  smaller  scales  grows.  This  causes  isosurfaces  in  the  coldest  regions  of  the  preheat 
zone  to  be  wrinkled  on  progressively  smaller  scales  (see  fig.  2).  In  this  context,  convergence  of  the  full  system  to  a 
unique  solution  in  the  purely  numerical  sense  would  not  be  expected  in  such  turbulent  flows.  At  the  same  time,  the 
observed  convergence  of  all  key  characteristics  of  the  turbulent  flame  simply  means  that  such  small  scale  motions  in 
cold  fuel  do  not  affect  the  evolution  of  the  flame  brush. 

Based  on  these  observations,  we  conclude  that  a  minimum  resolution  of  16  cells  per  6l  is  required  to  capture  the 
evolution  of  the  turbulent  flame  adequately.  Lower  resolutions  tend  to  exaggerate  both  the  values  and  the  degree  of 
variability  of  all  major  quantities,  primarily  the  turbulent  flame  width  and  speed. 

A  minimum  resolution  of  16  cells  is  substantially  higher  than  the  resolution  of  4  cells  per  6l  that  is  required  to 
reproduce  with  high  accuracy  the  behavior  of  the  planar  laminar  flame  [19].  This  suggests  that  the  need  for  higher 
resolution  arises  due  to  the  multidimensional  effects  and,  in  particular,  due  to  the  highly  curved  geometry  of  individual 
flame  sheets.  In  the  presence  of  high-speed  turbulence,  the  flame  is  often  folded  with  a  curvature  radius  comparable 
to  the  laminar  flame  width  (see  fig.  2).  In  such  cusps,  the  Cartesian  mesh  invariably  introduces  some  degree  of 
unphysical  anisotropy  to  the  thermal  flux.  This  effect  is  not  present  in  the  case  of  a  planar  laminar  flame.  In  the 
turbulent  flow,  however,  it  introduces  errors  which  at  low  resolutions  can  become  significant  in  comparison  with  the 
overall  solution  errors  that  occur  regardless  of  the  flame  geometry.  Consequently,  higher  resolution  is  necessary  to 
minimize  this  effect.  In  order  to  relax  this  minimum-resolution  threshold  and  bring  it  closer  in  line  with  that  required 
for  a  planar  laminar  flame,  it  is  necessary  to  improve  the  method  currently  implemented  in  Athena-RFX  [19]  in  terms 
of  the  multidimensional  coupling  of  diffusive  fluxes  as  well  as  their  coupling  with  the  hydrodynamic  fluxes.  This  is 
the  subject  for  further  work. 

4.4.  Implications  for  an  actual  stoichiometric  H2-air  mixture 

The  physical  model,  used  in  this  work,  relies  on  the  first-order  Arrhenius  kinetics  and  a  simplified  reaction- 
diffusion  model  of  [17]  for  stoichiometric  H2-air  mixture.  The  low  computational  cost  of  such  approach  makes  the 
large-scale  multidimensional  simulations,  such  as  the  ones  presented  in  this  work,  computationally  feasible.  More 
importantly,  however,  it  also  allows  us  to  accurately  reproduce  the  key  characteristics  both  of  laminar  flames  and 
multidimensional  detonations  in  stoichiometric  H2-air,  such  as  the  laminar  flame  width,  speed,  detonation  velocity 
and  detonation  cell  size,  as  well  as  the  dependence  of  these  quantities  on  pressure  and  temperature. 

In  order  to  investigate  the  role  of  small-scale  turbulence  on  the  flame,  we  did  not  include  physical  viscosity  in  our 
model  and  instead  relied  on  numerical  viscosity  to  provide  kinetic-energy  dissipation.  As  a  result  of  this,  however, 
our  simulations,  strictly  speaking,  do  not  describe  the  behavior  of  an  actual  stoichiometric  H2-air  mixture. 

It  was  discussed  in  §  2.4  that  fully  resolved  DNS  models  of  the  interaction  of  high-speed  turbulence  with  flames 
are  presently  unfeasible  because  the  Kolmogorov  scale  in  such  regimes  is  much  smaller  than  the  laminar  flame  width. 
In  light  of  this,  what  can  our  results  say  about  the  high-speed  turbulent  combustion  in  the  presence  of  the  actual 
H2-air? 

We  have  shown  that  scales  T  (5l  do  not  play  pronounced  role  in  determining  the  dynamics  and  properties  of  the 
turbulent  flame  brush.  This  suggests  that  turbulent  motions  on  yet  even  smaller  scales,  which  would  be  present  in  the 
cold  H2-air  mixture,  should  not  be  able  to  change  the  results  either  qualitatively  or  quantitatively. 

At  the  same  time,  the  evolution  of  the  turbulent  flame  is  completely  determined  by  scales  A  >  Si.  In  all  our 
simulations,  the  numerical  Taylor  scale  is  smaller  than  the  actual  Taylor  scale  in  the  product,  and  it  is  smaller  than 
the  actual  Taylor  scale  in  fuel  in  the  highest  resolution  case  S3.  Therefore,  motions  on  all  large  scales  are  properly 
captured  in  our  calculations.  This  fact,  along  with  the  negligible  role  of  small  scales,  suggests  that  our  models  can  be 
viewed  as  representative  of  the  actual  behavior  of  stoichiometric  H2-air. 
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The  only  exception  concerns  the  wrinkling  scale  of  the  flame  surface  on  the  fuel  side.  The  presence  of  turbulent 
motions  on  much  smaller  scales  would  result  in  substantially  finer  wrinkling  than  shown  in  fig.  2  even  in  our  highest- 
resolution  simulation.  At  the  same  time,  the  flame  surface  on  the  product  side  is  representative  of  that  in  the  actual 
reactive  mixture. 

Based  on  these  considerations,  our  results  suggest  that  the  classification  of  combustion  regimes  of  stoichiometric 
H  2-air  on  the  basis  of  the  internal  flamelet  structure  and,  in  particular,  the  degree  of  the  reaction-zone  broadening  does 
not  accurately  reflect  the  actual  process  of  turbulence-flame  interaction.  In  particular,  substantial  flame  broadening 
by  turbulence  appears  to  be  unlikely  for  hydrogen  in  all  subsonic  regimes.  In  our  simulations,  the  turbulent  Mach 
number  in  cold  fuel  is  Mqf  -  0.25.  Therefore,  in  order  for  the  turbulence-flame  interaction  to  remain  in  the  subsonic 
regime,  the  turbulent  velocity  can  be  increased  only  moderately.  This  does  not  appear  to  be  sufficient,  as  our  results, 
as  well  as  the  results  of  other  studies  [10,  15],  show  that  substantially  higher  turbulent  intensities  would  be  required 
to  disrupt  the  internal  flame  structure. 

While  faster  turbulence  would  be  more  likely  to  affect  the  internal  flame  structure,  it  would  also  cause  compress¬ 
ibility  effects  and  turbulent  heating  of  the  fuel  to  become  much  more  pronounced.  Resulting  higher  pressures  and 
temperatures  would  increase  the  laminar  flame  speed  while  decreasing  its  width,  thereby  reducing  the  broadening 
effect  of  turbulence.  It  is  also  not  clear  what  the  effects  of  the  stronger  compressive  component  of  the  velocity  field  in 
such  fast  turbulence  would  be.  In  particular,  local  compressions  and  rarefactions  can  also  modify  or  disrupt  the  local 
flame  structure.  Such  process  would  be  completely  distinct  from  the  traditionally  considered  action  of  the  vortical 
motions  associated  with  the  purely  solenoidal  part  of  the  turbulent  field.  We  consider  faster  turbulent  regimes  in  future 
work. 

Certainly,  the  simplified  reaction-diffusion  model  used  in  this  work  cannot  capture  the  full  complexity  of  hydrogen 
combustion.  While  we  do  not  expect  the  detailed  chemistry  to  change  presented  results  qualitatively,  it  is  nevertheless 
important  to  investigate  the  potential  effects  of  the  realistic  chemical  model.  This  is  the  subject  for  future  work. 

We  conclude  with  the  question  about  the  universality  of  the  observed  behavior.  Is  it  possible  to  identify  a  set  of 
parameters  which,  for  the  given  reactive  mixture  and  flow  conditions,  would  uniquely  specify  the  degree  of  flame 
broadening  and  the  role  of  small-scale  turbulence?  Are  the  given  ratios  of  the  turbulent  intensity  and  the  system  size 
to  the  laminar  flame  speed  and  width,  which  were  used  in  our  simulations,  sufficient  to  apply  the  obtained  results  to 
other  fuels? 

Our  results  appear  to  be  in  general  agreement  both  with  experimental  and  numerical  studies  [10,  15]  ofhigh-speed 
turbulent  combustion  in  such  different  systems  as  the  flat  thermonuclear  flame  in  degenerate  matter  and  the  jet-burner 
flame  in  compressed  natural  gas.  In  particular,  it  appears  that,  for  a  broad  class  of  reactive  mixtures,  the  effects  of 
turbulence  on  the  internal  flame  structure  are  significantly  suppressed,  and  disrupting  the  flame  requires  substantially 
higher  turbulent  intensities  than  can  be  predicted  based  on  simple  analysis. 
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A.  Isosurface-based  method  for  flamelet  structure  reconstruction 

Here  we  consider  in  detail  properties  of  the  method  that  we  use  to  reconstruct  the  average  flamelet  structure  in  the 
turbulent  flame  brush,  as  was  discussed  in  §  3.2. 

First  consider  a  planar  laminar  flame  normal  to  z-axis.  In  this  flame  all  isosurfaces  Si  are  exactly  parallel.  A,  = 
Const,  and  eq.  (13)  reduces  to  rji+i  =  z(T,+i)  -  z(T,)  where  z(Y)  is  based  on  the  exact  laminar  flame  profile  Y{z). 
Therefore,  Y(ri)  =  Y(z).  Hereafter,  for  brevity  we  only  consider  Y  and  the  same  reasoning  applies  to  T.  Imagine  now 
that  we  deform  that  planar  flame  without  stretching  it,  and  we  ensure  that  at  each  point  curvature  radius  is  much  larger 
than  the  flame  width.  In  this  case,  all  isosurfaces  remain  exactly  parallel.  A,  remains  constant,  and  again,  regardless 
of  how  complex  the  deformation  is,  this  method  will  recover  the  exact  laminar  flame  profile  Yirf)  —  Y(z). 

Next  consider  a  real  flame  subject  to  the  action  of  turbulence.  In  this  case  there  is  no  reason  to  assume  a  priori 
either  that  isosurfaces  Si  and  Si+i  are  parallel  or  that  their  surface  areas  are  equal.  Let  us  first  define  the  distance 
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d(p,  Si+\)  between  a  point  p,  located  on  isosurface  Si,  and  isosurface  >Si+i  as 

d{p,Si+\}^  min  Wp-p'h,  (14) 

where  ||.||2  is  the  Euclidean  norm.  Then  we  define  the  maximum  distance  between  the  isosurfaces  Si  and  >S,+i  in 
a  traditional  way  as  the  Hausdorff  distance 

^max  =  dmaxiSi,  Si+i)  =  max  |  max  (c/(p,  ^max  *Si))|.  (15) 

Similarly  we  define  the  minimum  distance  ^min  between  Si  and  >S,+i 

^min  =  dminiSi,  Si+i)  =  min  I  mm  (d(p,  *Si+i)),  ^mm  (d(p',  *Si))|.  (16) 

Using  these  definitions,  the  following  proposition  can  be  made: 


As  n  ^  oQ  and  AT  s  —  T,-  — >  0,  then  »S,+i  — >  Si  in  the  sense  that  as  ^min  0,  ^max  ^min- 

A  corollary  of  the  above  proposition  also  is  the  statement  that  A,  — >  A,+i  as  T,  — >  T,+i . 

In  order  to  prove  this  proposition  we  first  make  several  observations.  Isosurfaces  cannot  intersect  each  other,  since 
that  would  lead  to  infinite  gradients  and  unphysical  conditions.  Second,  since  the  evolution  of  Y  is  governed  partly  by 
diffusion  and  conduction,  and  there  are  no  shocks  in  the  system,  the  distribution  of  Y  is  smooth  and  continuous.  More¬ 
over,  Y  cannot  be  exactly  constant  over  any  extended  region  since  advective  and  diffusive  processes  in  an  unsteady 
turbulent  flow  would  inevitably  create  spatial  variations  in  the  distribution  of  Y  on  all  scales.  Therefore,  isosurfaces 
corresponding  to  all  values  Yi  <  Y  <  T,+i  are  all  bounded  by  the  isosurfaces  Si  and  >S,+i ,  they  are  smoothly  distributed 
throughout  the  volume  bounded  by  Si  and  >S,+i,  and  as  T,  — »  T,+i,  — >  0. 

Let  us  first  prove  that  a  value  Y'  e  (T,,  T,+i)  can  always  be  found  such  that  for  its  isosurface  S'  the  following  is 
true:  <  ^^ax  -  ^Lx-  Here  =  d^axiS' ,Si^\)  and  =  c/„,„(*S',*S,+i)  by  analogy  with  eqs.  (15)-(16). 

Assume  the  contrary,  namely  that  for  all  values  of  Y'  from  the  interval  (T,,  T,+i)  we  have  >  ^max-^'max-  This 

means  that  as  Y'  — >  T,+i,  — >  0  and  >  ^^ax  -  ^min-  Therefore,  S'  isosurfaces  have  some  limiting  isosurface 

S*  with  the  same  value  of  Y*  =  T,+i  but  which  is  distinct  from  >S,+i  since  >  ^^ax  -  ^min-  Since  isosurfaces  do  not 
intersect,  this  means  that  no  other  isosurfaces  pass  through  the  volume  bounded  by  S*  and  >S,+i.  Consequently,  S* 
and  >S,+i  bound  a  region  of  constant  Y.  However,  we  discussed  above  that  such  regions  cannot  exist.  Therefore,  we 
arrive  at  a  contradiction  and  our  assumption  was  wrong,  which  proves  our  initial  statement. 

It  follows  then,  that  there  exists  a  value  Y"  e  (T',  T,+i)  such  that  now  for  its  isosurface  S"  we  have  £'  ■  -  £"  ■  < 
^'max  ~  ^max-  HY  induction  this  proves  the  original  proposition.  A  similar  result  can  be  proven  for  temperature. 

In  practice,  this  result  means  the  following.  The  value  of  ^„jax  -  ^min  is  the  measure  of  the  variation  of  Si  with 
respect  to  >S,+i .  As  we  select  finer  discretization  intervals  of  Y,  this  variation  decreases  while  ^min  — >  0.  Therefore,  the 
isosurfaces  become  more  and  more  parallel  to  each  other.  An  illustration  of  this  can  be  seen  in  fig.  4.  While  there  is 
very  little  in  common  between  the  Y  -  0.05  and  Y  =  0.95  isosurfaces,  Y  -  0.6  tends  to  follow  the  Y  =  0.05  isosurface 
much  more  closely. 

An  important  question  concerns  the  choice  of  the  number  of  discretization  intervals  n  in  practical  applications.  It 
follows  from  the  above  discussion,  that  in  order  to  maximize  the  accuracy  of  the  method,  n  must  be  chosen  as  large 
as  possible.  This  ensures  that  consecutive  isosurfaces  are  close  to  each  other  and  ^^ax  -  ^min  is  minimal,  or  ideally 
zero.  In  the  computational  domain,  the  minimum  spatial  scale  is  set  by  the  cell  size  Ax.  This  determines  the  minimum 
practical  separation  of  isosurfaces  Isosurfaces  with  smaller  separations  would  pass  through  the  same  cell,  which 
would  result  in  a  substantial  drop  in  accuracy  of  the  overall  method  since  all  flow  variables  are  piecewise  constant 
within  a  cell.  At  the  same  time,  if  two  isosurfaces  on  average  have  separation  of  a  few  Ax,  this  ensures  that  their 
surface  areas  are  close  in  value,  and  there  cannot  exist  an  intermediate  isosurface  which  is  substantially  more  or  less 
tightly  folded.  These  considerations  allow  us  to  determine  the  maximum  value  of  n. 

If  we  assume  that  the  maximum  gradient  in  temperature  in  the  domain  is  close  to  that  in  the  laminar  flame  profile 
{dT I dx^L^ffiax^  then 


Ar„ 


{dT  ldx)L,„ 


Ax. 


(17) 


Recalling  that  6l  —  (Tp  -  To)l(dT ldx)L,max,  we  can  rewrite  eq.  (17)  as 


n 


Tp-To 

6l 


Ax, 


(18) 
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where  Tmin  and  T^ax  are  the  lower  and  upper  bounds  of  the  discretized  temperature  range.  Since  Tmin  ~  Tq  and 
Tmax  -  Tp,we  find  that 

nmax  ~  (19) 

Ax 

A  similar  argument  applies  to  the  fuel  mass  fraction.  Therefore,  n  should  not  be  larger  than  the  number  of  grid  cells 
within  the  laminar  flame  thermal  width.  Typically,  we  find  the  choice  of  n  based  on  eq.  (19)  adequate  in  practical 
applications. 

We  found  that  adaptivity  of  the  algorithm,  which  implements  the  method  described  here,  is  essential  to  guarantee 
high  accuracy.  The  estimate  given  by  eq.  (19)  assumed  that  the  largest  gradients  present  in  the  system  are  well 
approximated  by  (dYldx)L,max  and  (dT ldx)L,max  in  the  laminar  flame  profile.  This  may  not  be  the  case.  Therefore,  it  is 
important  for  the  algorithm  to  be  able  to  decrease  n  if  larger  gradients  are  encountered.  On  the  other  hand,  the  choice 
of  n  in  eq.  (19)  is  based  on  the  thermal  flame  width.  The  full  flame  width  is  larger  than  Sp-  Therefore,  while  eq.  (19) 
is  adequate  to  capture  flame  structure  in  the  steepest  region,  the  parts  of  the  profile  close  to  the  extreme  values  of  the 
discretized  range  may  be  under-resolved  (see  fig.  6).  In  those  regions  both  Y  and  T  vary  slowly  in  space.  Thus,  for 
consecutive  isosurfaces,  rji+i  »  Ax.  Ideally,  AT  and  AT  must  be  chosen  such  that  Ax  <  rji+i  <  aAx.  Typically,  we 
find  that  a  =  4  is  a  reasonable  choice. 

Then  the  adaptivity  of  the  algorithm  is  implemented  as  follows.  The  value  of  n  is  adjusted  to  ensure  that  both 
AT  and  AT  are  marginally  large  enough  to  prevent  situations  when  rji+i  <  Ax.  On  the  other  hand,  if  rji+i  >  aAx  is 
encountered,  AT  or  AT  for  that  interval  is  divided  by  two.  Subsequently,  rj  is  evaluated  separately  for  the  first  and 
second  half-intervals,  and  the  final  rji+i  is  the  sum  of  those  values.  If  any  of  the  intermediate  values  of  rj  is  larger  than 
aAx,  that  half-interval  is  again  subdivided  by  two  and  the  whole  procedure  proceeds  recursively  until  on  each  of  the 
substeps  the  condition  jj  <  aAx  is  satisfied.  This  ensures  that  t]  is  always  determined  for  isosurfaces  which  are  close 
to  each  other. 

Note,  that  this  algorithm  uses  the  laminar  flame  structure  only  to  provide  the  initial  guess  for  n.  Subsequently, 
it  adapts  the  discretization  intervals  for  T  and  T  only  based  on  the  actual  gradients  found  in  the  flow.  Therefore, 
in  principle,  the  information  about  the  laminar  flame  structure  is  not  required  at  all  and  eq.  (19)  is  intended  only  to 
facilitate  finding  the  correct  AT  and  AT.  At  the  same  time,  in  situations  when  the  internal  structure  of  the  turbulent 
flame  is  very  different  from  that  of  the  laminar  flame,  as  would  be  expected  in  the  case  of  distributed  burning,  eq.  (19) 
is  only  of  limited  benefit. 

A  key  limitation  of  this  method  is  that  this  procedure  can  be  applied  only  to  quantities  that  change  monotonically 
through  the  flame,  such  as  T,  T,  or  p.  In  particular,  this  method  cannot  be  used  to  directly  determine  the  distribution 
of  the  reaction  rate,  R.  Non-monotonic  behavior  of  a  quantity  means  that  there  is  an  inherent  degeneracy  in  its 
distribution  inside  the  flame.  Consider  a  planar  laminar  flame  (e.g.,  see  fig.  6).  Two  distinct  points  in  the  reaction 
rate  profile  have  the  same  value  of  R.  Therefore,  there  would  be  two  isosurfaces  Si  and  >S,+i  for  each  value  of  /?,  and 
Ri+\.  Volume  V,  would  be  the  total  volume  bounded  by  both  pairs  of  isosurfaces  and  area  A,  would  be  the  combined 
surface  area  of  both  isosurfaces  Si.  There  is  no  mechanism  in  this  method  to  distinguish  the  contributions  of  each  pair 
of  isosurfaces  both  into  A,  and  V,.  Consequently,  77,  found  using  eq.  (13)  would  have  no  meaning. 

In  general,  we  find  that  the  method  described  above  performs  best  for  quantities  that  predominantly  change  within 
the  volume  of  the  flame  brush  and  which  do  not  vary  at  all  outside  the  flame  brush,  or  they  vary  only  slightly  around 
the  limiting  values  of  the  discretized  range.  This  is  the  case  for  T  and  T .  The  accuracy  is  the  highest  for  T  which  is 
exactly  0  or  1  outside  the  flame  brush.  Therefore,  results  are  not  contaminated  by  the  isosurfaces  which  pass  outside 
the  burning  region.  For  temperature,  the  accuracy  is  marginally  lower  due  to  the  fluctuations  in  the  turbulent  field 
away  from  the  flame  brush,  which  can  cause  values  of  T  to  become  slightly  below  the  upper  bound  or  above  the  lower 
bound.  As  a  result,  isosurfaces  would  pass  through  such  regions  which  are  not  part  of  the  flame  brush.  We  find, 
however,  that  this  typically  introduces  small  errors  only  near  the  extreme  values  in  the  temperature  profile  (see  fig.  6). 
At  the  same  time,  for  a  quantity  like  density,  the  accuracy  drops  to  the  point  that,  while  we  were  able  to  obtain  the 
profiles,  we  find  that  they  cannot  be  reliably  used  to  analyze  the  density  structure  of  the  flamelets.  This  is  primarily 
due  to  the  fact  that  variations  in  density  outside  the  flame  brush  are  large  enough  to  be  significant  in  comparison  with 
the  change  in  p  inside  the  flame  brush  itself.  Therefore,  a  large  portion  of  the  profile  ends  up  being  contaminated  by 
the  contributions  from  the  regions  which  are  outside  the  flame  brush  but  which  happen  to  have  values  of  p  that  fall 
within  its  discretized  range. 

Finally,  this  method  does  not  provide  any  means  to  relate  the  distributions  of  the  obtained  quantities,  e.g.,  T  and 
T .  In  a  planar  laminar  flame  structure,  for  instance,  a  given  value  of  T  uniquely  defines  its  position  z  in  the  profile 
and  thus  uniquely  determines  the  value  of  T  =  T{z).  This  cannot  be  done  for  the  profiles  obtained  using  the  method 
described  above.  Recall  that  tj  is  not  a  spatial  coordinate  but  simply  the  distance  between  consecutive  isosurface 
values  and  it  is  assumed  that  rjo  -  0.  Thus,  values  of  tj  for  T  and  T  are  distinct  and  unrelated  and  other  arguments 
must  be  invoked  in  order  to  relate  their  distributions. 
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Table  1 .  Input  Model  Parameters  and  Computed  Laminar  Flame  Properties 


To 

293  K 

Initial  temperature 

Pq 

1.01  X  10®  erg/cm^ 

Initial  pressure 

Po 

8.73  X  10-“*  g/cm^ 

Initial  density 

r 

1.17 

Adiabatic  index 

M 

21  g/mol 

Molecular  weight 

A 

6.85  X  10*^  cmVg-s 

Pre-exponential  factor 

Q 

46.37  RTo 

Activation  energy 

q 

43.28  RTo  /  M 

Chemical  energy  release 

KO 

2.9  X  10^®  g/s-cm-K’^ 

Thermal  conduction  coefficient 

Do 

2.9  X  10-®  g/s-cm-K" 

Molecular  diffusion  coefficient 

n 

0.7 

Temperature  exponent 

Tp 

2135  K 

Post-flame  temperature 

Pp 

1.2  X  10^^  g/cm^ 

Post-flame  density 

Sl 

0.032  cm 

Laminar  flame  thermal  width 

Sl 

302  cm/s 

Laminar  flame  speed 

Table  2.  Parameters  of  Simulations^ 


Sl 

S2 

S3 

Description 

2) 

64  X  64  X  1024 

128  X  128  X  2048 

256  X  256  X  4096 

Domain  grid  size 

Da 

1  X  1  X  16 

Domain  aspect  ratio 

L 

0.259  cm  =  8di 

Domain  width,  energy-injection  scale 

Ax 

4.05  X  10“^  cm 

2.02  X  10“^  cm 

1.01  X  10-2  cm 

Cell  size 

Ax 

8 

16 

32 

di/Ax 

Zt.O 

1.95  cm  =  7.52L 

Initial  flame  position  along  z-axis 

£ 

1.26  X  10®  erg/cm^-s 

Energy-injection  rate 

Us 

4.53  X  10^  cm/s  =  15Sl 

Turbulent  velocity  at  scale  Sp 

U 

9.07  X  10^  cm/s  =  305^ 

Turbulent  velocity  at  scale  L 

^  rms 

1.04  X  10®  cm/s  =  34.48Sl 

Turbulent  r.m.s.  velocity 

u, 

5.60  X  10^  cm/s  =  18.54Sl 

Integral  velocity 

i 

6.04  X  10-2  cm=  1.87di 

Integral  scale 

'^ed 

2.86  X  10-®  s 

Eddy  turnover  time,  L/  U 

Ugn 

3.0T,d 

3.0T,d 

2.0T,d 

Ignition  time 

hotal 

16.0Ted 

Total  simulation  time 

La 

9.59  X  10-®  cm  =  2.96  x  IQ-^dz. 

Gibson  scale 

Da 

0.10 

Damkbhler  number,  {IIUi)I{6lIS p) 

ISdcif 

0.25 

Mach  number  in  fuel,  UICs,f 

Map 

0.09 

Mach  number  in  product,  UICs,p 

“  Parameters  common  to  all  simulations  are  shown  only  once  in  S2  column. 
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Table  3.  Time- Averaged  Width  and  Speed  of  the  Turbulent  Flame  Brush^ 


St  I5l 

0{StI6l) 

StISl 

O(StISt) 

SI 

16.13 

6.09 

S2 

14.86 

1.96 

4.50 

2.29 

S3 

14.42 

4.09 

“  Time-averaging  for  all  variables  is  performed 
over  the  time  interval  (2  -  16)Tej. 
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Figure  1  Instantaneous  energy  spectra  for  simulations  SI  (black),  S2  (red),  and  S3  (green)  at  a  time  immediately  prior 
to  ignition.  Shaded  regions  illustrate  scales  associated  with  thermal  width  Sl  and  full  width  26l  of  the  laminar  flame 
(cf.  fig.  6).  Vertical  dashed  lines  show  wavenumbers  corresponding  to  Taylor  microscales  in  product,  Ap,  and  fuel, 
Ap,  as  well  as  Kolmogorov  scale  in  product,  rjp,  based  on  the  value  of  the  viscosity  coefficient  discussed  in  §  2.4.  The 
Kolmogorov  scale  in  the  fuel,  rjp  -  1.18  x  10“^  cm,  is  located  outside  the  range  of  the  graph. 
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Figure  2  Flame  brush  structure  in  simulations  SI  (top  row),  S2  (middle  row),  and  S3  (bottom  row).  Shown  is  fuel 
mass  fraction  at  f  =  l^Ted-  Left  panels  show  view  from  the  product  side,  right  panels  show  view  from  the  fuel  side. 
Bounding  isosurfaces  represent  Y  -  0.05  and  Y  -  0.95.  The  thin  black  line,  corresponding  to  T  =  0.6,  marks  the 
boundary  between  the  preheat  and  reaction  zones.  The  thin  white  line,  corresponding  to  T  =  0.2,  shows  the  location 
of  the  peak  reaction  rate  (cf.  fig.  6). 
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Figure  3  (a)  Evolution  of  the  turbulent  flame  thickness  6t  normalized  by  5l-  Note,  that  domain  width  L  -  86l.  (b) 
Evolution  of  the  turbulent  flame  speed  S  j  normalized  by  5  £.  In  both  panels,  black  lines  correspond  to  simulation  S 1 , 
red  to  S2,  and  green  to  S3. 


Eigure  4  Isosurfaces  of  the  fuel  mass  fraction  in  simulation  S2  at  f  =  (cf.  fig.  2,  middle  row,  left  panel). 

Isosurface  values  are  0.05  (red),  0.6  (green),  0.95  (blue).  Red  and  green  isosurfaces  bound  the  flamelet  reaction  zone. 
Green  and  blue  isosurfaces  bound  the  preheat  zone.  The  zo.mm  and  zi^max  mark  the  flame-brush  bounds.  The  zo.max  and 
Zi,mm  indicate,  respectively,  the  maximum  extents  of  product  and  fuel  penetration  into  the  flame  brush. 
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Figure  5  Temperature  structure  of  the  flame  brush  in  simulation  S2  at  f  =  l3Ted  (cf.  fig.  2,  middle  row).  Left  panel 
shows  view  from  the  product  side,  right  panel  show  view  from  the  fuel  side.  Bounding  isosurfaces  represent  T  =  400 
K  and  T  -  2060  K.  Thin  black  line,  corresponding  to  T  =  1035/r,  marks  the  boundary  between  the  preheat  and 
reaction  zones,  while  thin  white  line,  corresponding  to  T  =  16^QK,  shows  the  location  of  the  peak  reaction  rate  (cf. 
fig.  6).  The  colormap  is  on  a  logarithmic  scale. 
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x/5j^ 


Figure  6  Time-averaged  flamelet  structure  in  the  turbulent  fl ame  brush.  Increasing  curves  represent  fuel  mass  fraction, 
decreasing  curves  represent  temperature.  Temperature  profiles  are  normalized  hy  Tp  (see  table  1).  Dashed  lines  are 
exact  solutions  for  the  planar  laminar  flame.  Black  lines  correspond  to  simulation  SI,  red  to  S2,  and  green  to  S3. 
Circles  represent  calculated  values  and  solid  lines  are  the  Akima  spline  fits.  Time  averaging  is  performed  over  the 
time  interval  (2  -  \6)Ted-  Open  blue  circles  represent  flamelet  structure  obtained  in  S3  using  half  the  time  averaging 
interval,  i.e.,  (9  -  \6)Ted-  Shaded  cyan  regions  show  the  range  of  variability  of  individual  profiles  in  S3  within  the 
time-averaging  interval.  Shaded  gray  region  shows  the  distribution  of  the  reaction  rate,  R,  in  the  exact  planar  laminar 
flame  solution  normalized  by  its  peak  value,  Rmax  -  9.5  x  10"^  s  ' . 
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a:  /  5^  (Laminar  flame  profiles) 

-1.5  -1.0  -0.5  0.0  0.5  1.0  1.5  2.0 


Figure  7  Instantaneous  average  distributions  of  the  fuel  mass  fraction,  Y,  temperature,  T,  and  the  reaction  rate,  R,  in 
the  simulation  S3  at  time  f  =  12. Each  profile  represents  the  distribution  of  a  given  quantity  along  the  z-axis, 
i.e.,  the  longest  dimension  of  the  domain,  with  each  point  being  the  average  over  the  x  -  y  plane.  For  comparison, 
corTesponding  profiles  for  the  planar  laminar  flame  are  superimposed  with  dashed  lines.  The  bottom  coordinate  scale 
is  for  the  turbulent  flame  while  the  top  scale  is  for  the  laminar  flame.  Note  the  difference  in  the  range  of  these  two 
scales.  In  particular,  the  turbulent  flame  profiles  are  more  than  a  factor  of  4  wider.  Similar  to  fig  6,  T  and  R  are 
normalized  by  their  respective  peak  values  in  the  exact  planar  laminar  flame  solution,  namely  Tp  (see  table  1)  and 
Rmax  -  9.5  X  10"^  s  '.  Positions  in  the  graph  are  shifted  so  that  the  coordinate  of  the  Y  -  0.5  point  in  the  fuel  mass 
fraction  profiles  is  0.0  as  in  fig.  6. 
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